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I.  INTRODUCTION 


In  the  first  report  [1]  on  discrete  convolution  method 
for  solving  some  large  moment  equations,  three  types  of  one 
and  two  dimensional  problems  were  solved.  In  this  report, 
we  solve  three  more  types  of  one  and  two  dimensional 
problems.  They  are 

(i)  scattering  from  a  helix 

(ii)  planar  arrays  with  antenna  elements  arranged  in 
triangular  pattern,  solved  using  one  expansion 
function  per  element 

(iii)  planar  arrays  with  antenna  elements  arranged  in 
triangular  pattern,  solved  using  three  expansion 
functions  per  element 


As  in  the  first  report,  the  computing  time  measurements 
(made  on  a  KL/lO  machine)  and  the  number  of  iterations 
needed  for  the  given  accuracy  are  listed.  In  addition,  the 
array  factor  of  the  planar  arrays  are  given. 


II.  SAMPLE  COMPUTATIONS  AND  COMMENTS 

The  "one"  dimensional  problem  is  scattering  from  a 
helix.  Fig.  1  shows  the  problem  of  scattering  from  a  helix. 
The  MOM  formulation  of  the  helix  problem  requires  that  the 
helix  be  subsectioned  into  equal  length  segments.  If  we 
number  the  helix  segments  so  that  the  segment  numbers  are  in 
consecutive  increasing  order  from  top  to  bottom,  then  it  is 
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Fig.  1. 


The  helical  wire  scatterer 
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Fig.  2. 


The  helical  wire  scatterer  with  gaps 
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easy  to  see  that  the  mutual  coupling  between  segments  as 
given  by  the  MOM  impedance  matrix  Z  is 

Z  =  Z (m-n)  (1) 
mn  ' 

i.e.,  the  value  is  dependent  only  on  the  difference  between 
segment  numbers.  Therefore  the  Z  matrix  will  be  Toeplitz  if 
the  helix  is  a  complete  helix  as  shown  in  Fig.  1.  The  Z 
matrix  will  be  non-Toeplitz  if  the  helix  has  gaps  in  between 
as  shown  in  Fig.  2.  It  is  apparent  from  the  discussions  of 
other  one  dimensional  problems  that  both  problems  can  be 
solved  using  the  one  dimensional  DCM  technique.  For  the 
helix  with  gaps  all  we  need  is  to  insert  phantom  elements  as 
shown  in  Fig.  2. 

Table  1  gives  the  number  of  iterations  needed  to  get 
the  required  degree  of  accuracy  for  the  helix  problem.  We 
can  see  that  the  number  of  iterations  needed  is  practically 
independent  of  the  length  of  the  helix. 

The  problem  of  a  planar  array  with  antennas  arranged  in 
triangular  patterns  instead  of  rectangular,  can  also  be 
formulated  as  a  two  dimensional  convolution  equation  by 
adding  phantom  elements  (as  shown  in  Fig.  4),  to  make  a 
parallelogram.  The  triangular  pattern  arrangement  is  shown 
in  Fig.  3. 

The  MOM  formulation  using  one  expansion  per  antenna 
then  gives  a  block  Toeplitz  matrix  which  can  be  solved  using 
two  dimensional  DCM.  However,  it  cannot  be  solved  using  the 
block  Toeplitz  method  since  the  field  on  the  phantom 
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Table  1.  Results  for  some  helical 
Radius  Pitch  Number  Number  of 
of  turns  segments 

.125  .30  6  120 

.125  .30  6  120 

.15  .50  6  120 

.15  .30  6  120 

.15  .30  8  160 

.15  .30  8  160 


scatterer  problems 


Number  of 

Field 

Last 

i tera tions 

error 

current 

change 

(%) 

(%) 

22 

.693 

.823 

.0498 

.1 

*• 

23 

1.01 

1.11 

.141 

.113 

8 

.812 

1.39 

.0504 

.20  2 

11 

.0562 

.0958 

.0036 

.0142 

10 

1.56 

1.17 

.177 

.248 

16 

.0668 

.0496 

.0076 

.0109 

12 

.568 

1.21 

.0375 

.0778 

16 

.0697 

.  153 

.0046 

.0097 

13 

.453 

.883 

.0572 

.0467 

# 

16 

.0934 

.184 

.0118 

.0097 
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Table  1  •  continued 

Radius  Pitch  Number  Number  of  Number  of  Field  Last 


of  turns  segments 

iterations  error  < 

current 

change 

(%) 

(%) 

* 

.15 

.30  8  160 

11 

1.07 

.785 

.104 

.141 

* 

16 

.0821 

.0602 

.0082 

.0114 

.15 

.30  16  320 

13 

.395 

.659 

.0154 

.0295 

16 

.0848 

.139 

.0033 

.00634 

.25 

.50  8  160 

11 

.537 

1.86 

.0282 

.2 

16 

.0495 

.  163 

.0026 

.0184 

Here , 

Radius  is  the  radius  of  the 

helix  in 

wavelengths 

Pitch  is  the  pitch  of  the  helix  in  wavelengths 
the  wire  radius  in  all  cases  is  .006  wavelengths 
first  entries  for  each  problem  in  the  last  two  columns 
are  maximum  values  and  second  entries  are  average  values 
t  indicates  that  polarity  of  impressed  E  field  is  not  as 
shown  in  Fig.  10  but  is  z  directed 

*  indicates  that  the  impressed  field  is  not  from  the 
direction  shown  in  Fig.  10  but  is  from  z  direction 


ig.  3.  The  triangular  pattern  arrangement. 


/Y 


s 


\ 


Phantom 

flfmemtS 


A  planar  array  with  antennas  arranged  in 
angular  pattern. 
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elements  are  unknown.  .  Therefore,  only  LU  decomposition  or 
two  dimensional  DCM  can  be  used.  For  a  large  array,  DCM 
will  be  considerably  faster.  .Since  the  current  on  each 
antenna  is  not  symmetric,  using  three  expansion  functions 
per  antenna  gives  a  much  more  accurate  result  and  needs  five 
times  more  computing  time  for  the  DCM.  With  LU 
decomposition  method  computing  time  will  go  up  twenty  seven 
times  the  already  large  value. 


Table  2  lists 
iterations  needed 
function  per  antenna 
time  and  number  o 
using  three  expansio 
the  problems  are 
antennas  one-quarter 
ground  plane.  The 
wavelength  in  either 


the  computing  time 
for  the  DCM  solution  usi 
element.  Table  3  lists 
f  iterations  needed  for 
n  functions  per  antenna 
for  planar  arrays  with 
wavelength  in  front  o 
seperation  between  ante 
d  i  rection . 


The  graphs  given  in  Figs.  6,  7,  8,  9,  1 
the  arrays  factors  for  the  planar  arrays 
pattern  arrangement  solved  by  the  DCM  using 
function  per  antenna  element.  The  arr 
computed  in  the  plane  perpendicularly  bisect 
shown  in  Fig.  5.  Angle  measurements  are  a 
the  figures,  the  solid  lines  give  the  array 
solutions  which  take  the  mutual  coupling 
into  account  and  the  dashed  lines  are  for 
solutions  which  do  not  take  the  mutual 


and  number  of 


ng  one 

expansion 

the 

computing 

the  DCM 

solution 

eleme 

nt.  All 

0.48  wavelength 

f  the 

infinite 

nnas  is 

one-half 

0 ,  and 
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are 

wi  th  t 

r iangular 

one 

expan 

sion 

ay  f  ac 

tors 

are 

ing  the 

arra 

y  as 

s  shown 

.  In 

all 

factors 

for 

the 

between 

antennas 

the 

ideal 

i  zed 

coupl 

ing 

into 
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Table 

2.  Results 

for 

some  planar 

arrays  with 

triangular 

pattern 

arrangement 

N 

Excitation 

I 

Computing 

Field  Error 

Current 

Time ( secs) 

(%) 

Change (%) 

12 

Uni fo rm 

5 

4 

.1789 

.664 

.0867 

.378 

Beam  steer 

6 

4 

.1124 

.595 

(45*  ) 

.0554 

.173 

76 

Uniform 

6 

20 

.1447 

.489 

.0283 

.116 

Beam  steer 

6 

20 

.2768 

.796 

(45*  ) 

.0333 

.098 

372 

Uniform 

5 

105 

.4950 

1.674 

.0467 

.174 

Beam  steer 

6 

105 

.3636 

.863 

(45*) 

.0163 

.0456 

Here, 


N  is  the  number  of 
I  is  the  number  of 
accuracy.  For  both 
change,  the  upper 


antennas  in  the  array 
iterations  needed  to  get 
field  error  and  (last)  c 
entry  is  the  maximum  and 


the  given 
ur rent 
the  lower 


entry  is  the  average 
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Table 

3.  Results 

for  some  planar 

arrays  with  triangular 

pattern 

arrangement  (Multiple 

Expansion  Solutions) 

N 

Excitation 

I  Computing 

Field  Error 

Current 

Time (secs) 

(%) 

Change ( %) 

76 

Uniform 

6  67 

.0401 

.3612 

.0027 

.0836 

Beam  steer 

6  67 

.078 

.6898 

(xz  45* ) 

.0032 

.0704 

Beam  steer 

6  67 

.0251 

.527 

(yz  135*  ) 

.0016 

.0493 

372 

Uniform 

5  165 

.1450 

1.298 

.0048 

.134 

Beam  steer 

5  165 

.2928 

2.1100 

(xz  45*  ) 

.0047 

.0994 

Beam  steer 

5  165 

.0706 

1.43 

(yz  135*) 

.0029 

.0748 

Here 

,  N  is  the 

number  of  antennas  in 

the  array 

I  is  the 

number  of  iterations  needed  to 

get  the  give 

accuracy . 

For  both  field  error 

and  (last)  current 

change,  the  upper  entry  is  the 

max imum 

and  the  lower 

entry  is 

the  average. 
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Fig.  5.  The  relative  position  of  the  plane  in  which 
the  array  factors  are  computed. 


ARRAY  FACTOR 
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o 
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13 


o 


Fig.  7  Array  factor  of  the  twelve  antenna  planar  array  with 
excitation  to  give  a  45  degrees  scan 
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o 
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16 


O 

o 


17 


o 

O 


Fig.  11  \rray  factor  of  the  three  hundred  and  seventy  two 
antenna  planar  array  with  excitation  to  give  a  45 
degrees  scan 
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account. 

As  we  can  see  from  the  graphs,  for  larger  arrays,  the 
main  beam  is  not  effected  by  the  mutual  coupling  but  the 
side  lobes  and  the  nulls  are  effected  strongly. 

The  graphs  given  in  Figs.  13,  14,  15,  15,  17,  18,  19, 
20,  21,  22,  23,  and  24  are  the  arrays  factors  for  the  planar 
arrays  with  triangular  pattern  arrangement  solved  by  the  DCM 
using  three  expansion  functions  per  antenna  element.  The 
array  factors  are  computed  in  the  xz  and  yz  planes 
perpendicularly  bisecting  the  array  as  shown  in  Fig.  12. 
Angle  measurements  are  as  shown.  In  all  the  figures,  the 
solid  lines  give  the  array  factors  for  the  solutions  which 
takes  the  mutual  coupling  between  antennas  into  account  and 
the  dashed  lines  are  for  the  idealized  solutions  which  do 
not  take  the  mutual  coupling  into  account. 

The  array  factor  is  defined  as  the  far  field  pattern 
divided  by  the  element  factor.  Therefore  for  the  DCM 
solution  using  one  expansion  function  per  element,  the 
element  factor  is  the  far  field  pattern  of  the  expansion 
function  and  so  the  array  factor  can  be  and  is  computed 
directly  from  the  solved  currents.  However,  with  three 
expansions  per  antenna  element,  the  current  distribution 
over  each  element  is  different  from  the  others.  Therefore, 
since  our  main  purpose  in  computing  array  factors  is  to 
compare  the  far  field  patterns,  we  define  a  normalized  array 
factor  as  the  total  far  field  pattern  divided  by  the  far 
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20 


o 

O 


Fi^.  '3  Army  factor  in  the  *?.  piano  of  tho  seventy  six 
antenna  planar  array  with  uniform  excitation 
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o 

O 


Fig.  \4  Arroy  factor  in  the  yz  plane  of  the  seventy  nix 
antenna  planar  array  with  uniform  excitation 
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o 
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Fiq.  15  Array  factor  in  the  xz  plane  of  the  seventy  six 

antenna  planar  array  with  excitation  to  give  a  45 
degrees  scan  in  the  xz  plane 
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o 
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24 


o 
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25 


o 

O 


FiQ.  18  Army  factor  in  the  yz  piano  of  tho  seventy  six 

antenna  planar  array  with  excitation  to  give  a  135 
decrees  scan  in  tho  yz  plane 
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o 

O 


Fig.  19  Array  factor  in  the  xz  plane  of  the  three  hundred 
and  seventy  two  antenna  planar  array  with  uniform 
excitation 
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CD 

o 


Fi<7-  20  Array  factor  in  the  yz  piano  of  the  three  hundred 
and  seventy  two  antenna  planar  array  with  uniform 
excitation 
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o 
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29 


o 

O 


Fiq.  12  Array  factor  in  the  yz  plane  of  the  throe  hundred 
and  seventy  two  antenna  planar  array  with 
excitation  to  <jive  4  5  decrees  scan  in  the  xz  plane 


30 


o 

O 


Fig.  23  Array  factor  in  the  xz  piano  of  tho  throe  hundred 
and  seventy  two  antenna  planar  array  with 
excitation  to  give  130  degrees  scan  in  the  yz  plane 
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field  pattern  of  an  element  with  only  one  expansionn 
function.  This  is  the  array  factor  used  in  the  graphs 
discussed  above. 

(pp)n 

(NAF)  n  =  -  (2) 

(EP)  x 

Here  (NAF)  is  the  normalized  array  factor  for  the  DCM 
solution  using  n  expansions  per  antenna  element.  (FP)  is 
the  far  field  pattern  for  the  DCM  solution  using  n  expansion 
functions  per  antenna  element.  (EP)^  is  the  element  pattern 
or  the  far  field  pattern  of  an  element  with  only  one 
expansion  function  per  element.  Since  both  (FP)n  and  (EP)^ 
goes  to  zero  when  the  angles  are  either  0  or  180  degrees, 
(NAF)  is  indeterminate  at  those  angles. 

As  we  can  see  from  the  graphs,  for  larger  arrays,  the 
main  beam  is  not  effected  by  the  mutual  coupling  but  the 
side  lobes  and  the  nulls  are  effected  strongly.  In  addition 
we  can  see  that  the  difference  between  the  actual  and  the 
ideal  array  factors  are  more  pronounced  for  the  45  degrees 
beam  steering  in  the  xz  plane.  This  is  expected  since  the 
coupling  between  the  antenna  elements  that  are  adjacent  in 
the  x-direction  is  stronger  than  the  coupling  between  the 
antenna  elements  that  are  adjacent  in  the  y-direction. 
Also,  we  can  see  that  the  array  factors  computed  from  the 
solutions  using  one  expansion  function  per  antenna  and  the 
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array  fractors  Computed  from  the  solutions  using  three 


expansion 

functi 

ons 

per 

antenna 

are 

surpr is  i 

ngly  close. 

This  ind 

icates 

that 

at 

least 

for  th 

e  given 

spacings  and 

antenna 

sizes. 

the 

sol 

utions 

using 

only  one  expansion 

function 

per  antenna 

may 

be  good 

enough 

for  most 

purposes . 

III.  CONCLUSIONS 


Three  more  types  o 
report  on  DCM  [1]  ar 
needed  are  found  to  be 
array  problems.  The  he 
but  the  number  of  itera 
accuracy  is  found  to 
length  of  the  helix. 


f  problems  sugge 
e  solved  and  the 
reasonably  smal 
lix  problems  requ 
tions  needed  for 
be  practically 


sted  in 
number  of 
1  for  th 
ired  more 
the  requir 
i ndepend 


the 
i  te 
e 

i  te 

ed 

ent 


first 

rations 

antenna 

rations 

degree  of 

of  the 


Three  other  major  types  of  problems  formulated  or 
suggested  in  Tl]  still  remains  to  be  solved.  They  are 

(i)  scattering  from  a  arbitrarily  shaped  planer 
conducting  surfaces  or  dielectric  sheets 

(ii)  scattering  from  a  imperfectly  conducting  or 
dielectric  bodies 

(iii)  planar  array  antenna  backed  by  a  finite  sized 
ground  plane 


The  major  question  here  is  wherether  or  not  convergence  will 
be  achieved  for  these  types  of  problems.  The  condition  for 
convergence  is  that  the  absolute  value  of  the  largest 
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eigenvalue  of  must  be 
are  some  practical  applic 
types  of  problems  are 
representative  problems  a 
extension  to  this  report. 


less  then  1  a 
ations  where 
useful.  Th 
re  suggested 


s  shown  in  [ 1 ] . 
the  solutions  o 
erefore,  solvi 


There 
f  these 
ng  some 


as  possible 


future 
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APPENDIX  A 

SPATIAL  DOMAIN  INTERPRETATION  OF  THE  SPECTRAL  ITERATION 

TECHNIQUES 


Both  the  spectral  solution  techniques (STD, SIT  etc.) 
[2],  [3],  [4],  [5]  and  the  Discrete  Convolution  Method  (DCM) 
[1]  use  the  discrete  Fourier  Transform  to  solve 
electromagnetic  fields  problems.  The  following  spatial 
domain  interpratation  of  the  spectral  solution  techniques 
may  give  a  better  understanding  of  the  differences  between 
the  methods. 

The  electric  field  integral  equation  for  a  perfectly 
conducting  scatterer  or  radiator  is  [2] 

(  G  *  J  )  =  r€S  (Al) 


Here  G  is  the  Green's  Dyadic,  J  is  the  surface  current 
density  on  surface  S  and  E.  is  the  tangential  component  of 
the  impressed  field.  We  can  extend  the  field  equation  over 
all  space  so  that  (Al)  becomes 


-*  i  — • 

-E  +  F 


( A2 ) 
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where  F  is  the  field  outside  the  surface  S.  The  Fourier 
tranformed  version  of  (A2)  reads 

/\>  ^  ^  ^ 

G(-O-)  j  (XL)  =  (.O.)  +  F(-n-)  (A3 ) 

-£_l-  =  (  toX,  co^,  uoz)  in  general 

A  formal  solution  to  (A3)  is 

J(Jt)  =  G_1(-0.)  l(-n.)  ( A4 ) 

where  E  (-&-  )  is  -E^  (-O. )  +F  (_QJ  . 

The  spectral  solution  techniques  essentially  consists 
of  iterative  techniques  that  find  F  and  thus  jf.  However 
since  f  cannot  be  found  in  closed  form,  the  solution  is  a 
numerical  solution.  Therefore  (A4)  and  hence  (A3)  is  solved 
only  for  a  finite  number  of  specific  discrete  frequencies. 
The  solution  we  get  is  an  exact  solution  of  (A3)  and  hence 
of  (Al)  only  if  ?  is  known  exactly  and  if  the  solution  is 
for  all  frequencies.  In  the  following  discussion  we  will 
set  aside  the  question  of  inaccuracies  due  to  inexact 
knowledge  of  F  because  it  is  a  question  of  convergence. 
Here  we  are  discussing  spatial  domain  interpretation  rather 
than  the  convergence  of  the  solution.  Therefore 

(a)  since  (A3)  is  satisfied  only  for  discrete  frequencies, 
we  are  sampling  in  frequency  domain  and  we  are  solving  the 
following  set  of  equations 
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G  (-n-)J  (-n.)  &  (Sl-SL^)  =  E(-0.)S  (XL-^1  )  ( A5 ) 


r l 

=  (60X, 

coY 

coz  ) 

n 

a  ' 

b ' 

c 

coX 

a 

,  X 

=  60  , 

0  ' 

•  •  • 

co  y 
b 

■ 

“i' 

•  •  • 

c 

“V 

•  •  • 

instead  of  (A3) . 

(b)  Since  we  cannot  solve  for  infinite  number  of  discrete 
frequencies,  we  solve  for  only  a  finite  number.  This 
amounts  to  truncating  in  frequency  domain.  If  we  denote  the 
truncating  function  by  H(-Tt),  we  are  solving  the  following 
set  of  equations 


H(-0-)G  (-H-  )  Ji-TL)  S  =  H  (XL)  E  (XL)  (A6) 

instead  of  either  (A3)  or  (A5). 

Before  we  give  the  spatial  domain  interpretation  of 
( A6 ) ,  we  will  first  examine  the  errors  due  to  above 
approximations  from  a  frequency  domain  perspective.  It  is 
well  known  [6]  that  sampling  in  the  spatial  domain  will  lead 
to  aliasing  in  frequency  domain  if  the  function  sampled  is 
not  limited  in  frequency.  By  reason  of  symmetry,  we  can 
easily  show  that  sampling  in  frequency  domain  will  lead  to 
aliasing  in  the  spatial  domain  if  the  function  sampled  is 
not  limited  in  space.  In  solving  (A3),  we  used  sampling  on 
both  sides  of  the  equation.  The  surface  current  J  is  indeed 
limited  in  space  (being  confined  to  conductor  surface  S)  but 
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the  £ield  E  is  not.  Therefore,  there  will  be  aliasing 
irrespective  of  the  sampling  rate.  However,  since  the  field 
E  decays  outside  the  surface,  aliasing  effect  can  be  reduced 
by  a  high  enough  sampling  rate  so  that  the  overlap  occurs  in 

-r 

the  region  where  E  has  already  decayed  to  a  small  value. 
Truncation  of  the  frequency  range  by  (the  windowing 
function)  on  both  sides  of  the  equation  leads  to  Gibb's 
phenomenon.  This  windowing  error  can  be  reduced  by  using 
different  windowing  strategies  discussed  in  [6],  A  good 
example  is  the  Hamming  windowing  function. 


We  will  now  look  at  the  above  errors  from  a  spatial 
domain  perspective.  In  spectral  solution  techniques, 
solution  to  ( A  6 )  is  found  by  using  the  Discrete  Fourier 
Transform.  This  is  possible  if-TL^'s  are  chosen  at  equally 
spaced  intervals.  It  can  be  proved  [6]  that  (A 6)  in  this 
case  is  equivalent  to  the  Discrete  Convolution  equation 


*  J (K)  =  E  (K) 

k  =0,1, 
x 

ky  =  0,1, 

k  =0,1, 
z 


I<  =  ( k  ,  k  ,  k  ) 
x  y  z 

,L 

,N 


(A7) 


Up  to  this  point,  we  have  been  discussing  the  general 
multi-dimensional  case.  From  here  on  however,  we  will 
confine  our  discussion  to  the  one  dimensional  case. 
Discussion  for  two  or  three  dimensional  problems  would  be 
similar.  (We  will  use  lower  case  indices  instead  of  higher 
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case,  to  indicate  one  dimensionality.)  Since  Ej  is  known 
exactly  and  F,  once  we  achieve  convergence,  is  also  known 
exactly,  E(k)  represents  the  sampled  value  at  kth  point  of 
the  exact  field  E  [2],  [51.  Also,  J  is  the  required 

solution  and  hence  J  (k)  represents  the  sampled  value  at  kfc^ 
point  of  the  current  J.  But  ^-(k)  is  not  the  sampled  value 
of  G.  Instead,  C^,  (k)  is  the  inverse  Discrete  Fourier 
Transform  (based  on  N+l  points)  of  G  ( to )  . 


Y,  G  (n)  e3  (  W4i  >  kn 


(A8) 


Therefore  Lj.(k)  is  the  inverse  Fourier  Transform  of  the 
sampled  and  truncated  (or  in  general  windowed)  G(tO). 

It  is  well  known  [6]  that  sampling  in  spatial  domain 
will  lead  to  a  function  in  frequency  domain  which  consists 
of  periodic  repetitions  of  the  original  frequency  domain 
function.  By  reason  of  symmetry  between  the  domains,  we  can 
easily  show  that  sampling  in  frequency  domain  will  lead  to  a 
function  in  spatial  domain  which  consists  of  periodic 
repetitions  of  the  original  spatial  domain  function. 
Therefore,  sampling  of  G  ( U> )  will  cause  the  inverse  of  the 
sampled  function  to  consist  of  periodic  repetitions  of 
the  original  spatial  domain  function  G.  This  will  cause 
overlaps  as  shown  in  Fig.  25  for  the  one  dimensional  case. 
But  G  is  the  Green's  Dyadic,  the  field  due  to  a  unit  impulse 
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Fig.  25  Relationship  between 

the  one  dimensional  case 


and  G  for 
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current  at  the  origin.  Therefore  C^.  g  is  the  overlapped 
approximate  field  of  a  unit  impulse  current  at  the  origin. 
Using  higher  sampling  rates  will  cause  the  rate  of 
repetition  to  decrease,  i.e.,  increase  X.  Since  we  are 
interested  in  solving  (A7)  for  only  a  finite  region  in 
space,  error  due  to  overlap  will  be  insignificant  if  the 
sampling  rate  is  high  enough. 


However,  in  addition  to  sampling,  we  also  truncate  the 
frequency  function.  This  truncation  will,  in  general,  be 

ry 

done  by  using  a  windowing  function  H  ( to  )  .  Therefore  G(k) 
can  be  thought  of  as  the  inverse  Fourier  Transform  of  the 
sampled  version  of  the  function  G(to)H(co).  G(co)  is  the 

Fourier  Transform  of  the  Green's  Dyadic  G,  the  field  due  to 

—  ^ 

a  unit  current  impulse  at  the  origin.  Therefore  G(co)h(<~j) 
is  the  Fourier  Transform  of  the  field  due  to  a  current 
distribution  H,  where  H  is  the  inverse  Fourier  Transform  of 


the  windowing  function  H(^).  Therefore  Lp(k)  is  the 
sampled  value  at  kfc^  point  of  the  field  due  to  the  current 
distribution  H.  Therefore,  the  result  of  the  convolution 


product  ^-*J  is  the  (overlapped)  approximate  field  due  to 

• — ^ 

the  sampled  values  of  the  current  J  multiplied  by  the 
current  distribution  H.  Since  a  convolution  product 
produces  a  finite  set  of  values,  each  one  corresponding  to  a 
point  in  space,  each  value  can  be  interpreted  as  the 
approximate  field  at  the  corresponding  point  in  space  due  to 
the  current  J,  expanded  in  terms  of  the  basis  or  expansion 
function  H.  Since  solving  (6)  amounts  to  matching  the  left 
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hand  side  and  the  right  hand  side  of  (A7)  [6]  ,  the  solution 
of  (A6)  can  now  be  seen  as  point  matching  the  field  E  with 
the  approximate  field  due  to  J  expanded  in  terms  of  the 
basis  or  expansion  function  "h.  If  the  error  due  to  overlap 
is  not  significant,  then  the  solution  of  (A6)  is  equivalent 
to  point  matching  the  known  field  E  with  the  field  due  to  J 
expanded  in  terms  of  the  basis  or  expansion  function  H. 
Therefore,  the  spectral  iteration  method  is  equivalent  to 
the  Method  of  Moments  [7] ,  if  we  use  the  inverse  Fourier 
Transform  of  the  windowing  function  H(co)  as  expansion 
function  and  point  match.  The  only  difference  is  that  with 
the  Method  of  Moments,  Gaussian  elimination  is  normally  used 
instead  of  iterative  techniques. 

The  Discrete  Convolution  Method  [1]  on  the  other  hand, 
is  the  iterative  solution  of  the  matrix  equation  formulated 
by  the  Method  of  Moments.  Therefore  the  difference  between 
the  Discrete  Convolution  Method  and  the  spectral  iteration 
techniques  is  that  with  the  former  method,  the  expansion 
function  is  chosen  explicitly  whereas  with  the  latter,  the 
expansion  function  is  chosen  implicitly;  i.e.  the  inverse 
Fourier  Transform  of  the  windowing  function. 
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APPENDIX  B 
COMPUTER  PROGRAMS 


The  computer  programs  given  in  this  section  are  written 
to  solve  the  three  additional  types  of  problems.  Although 
they  are  not  optimized  in  any  sense  of  the  word,  they  are 
wr itten  to  avoid  obvious  wastes  of  computing  time  and  memory 
space  and  is  fairly  efficient. 


I.  MAIN  PROGRAM  SEGMENT  AND  SUBROUTINES  FOR  THE  HELIX 
PROBLEM 

Both  the  main  program  segment  and  the  subroutine  CALZ 
are  from  [8] .  The  slight  modifications  in  the  main  program 
are  self-explaining  so  no  attempt  will  be  made  here  to  give 
a  description.  The  subroutine  HELIX  on  the  other  hand,  is 
written  to  compute  the  co-ordinates  of  the  helix  points  from 
the  given  radius,  pitch,  number  of  turns,  and  the  number  of 
points  for  which  the  co-ordinates  must  be  found. 


400 

c 
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500 

c 

THIS  IS  THE  MAIN  PROGRAM 

600 

c 

700 

COMPLEX  Z(3600)  ,0(400)  ,C(400)  ,E(3)  ,EI<2)  ,OV  (4) 

800 

COMPLEX  01,ZL{30) ,ZIN,TIN,V,ZC 

900 

COMPLEX  TOE{1024),EX(400) ,CUR  (1024) ,GOESS 

000 

COMMON  /G/GOESS,COR 

100 

COMMON  XX (800)  ,XI(800)  ,XZ(800)  ,TX(800)  ,TY(800)  ,TZ(800)  ,AL(800) 

200 

COMMON  T  (1600)  ,TP(1600)  ,BK,RAD2(21),L(11),LR  (21)  /COA/C 

300 

DIMENSION  PX  (900)  ,PY  (900)  ,PZ(900)  ,LL(11)  ,RAD(21)  ,IFP(60) 

400 

DIMENSION  LP  (30) 

500 

PI=3. 141592654 

550 

OPEN  (0NIT=21 ,FILE=* ARDATA.DAT’) 

600 

1 

READ  (1,101)  NW, NP,NR,BK,IFLAG 

700 

WRITE  (3,  102)  NW,NP,NR,  BK 

800 

LL  (NW+ 1) =801 

900 

LR  (NR+ 1)  =80 1 

000 

IF(IFLAG.NE.O)  GO  TO  300 

100 

READ  ( 1 ,  1  30)  (PX  (I)  ,I=1,NP) 

200 

READ(1,130)  (PY  (I)  ,1=  1 ,  NP) 

300 

R EAD  ( 1 ,  130)  (PZ  (I)  ,  1=1 ,  NP) 

400 

GO  TO  301 

500 

300 

CONTINOE 

510 

IF (IFLAG. NE. 1)  GO  TO  310 

520 

CALL  HELIX (PX,PY,PZ,NP) 

530 

GO  TO  301 

540 

310 

CONTINOE 

600 

PE AD  (1, 130) WIREL, SEGL 

700 

DO  501  1=1, NP 

800 

PX  (I) =WIREL 

900 

WIREL=WIREL+SEGL 

000 

PY (I) =0. 0 

100 

PZ  (I)  =0.0 

200 

501 

CONTINUE 

300 

301 

CONTINUE 

100 

WRITE  (3,  104)  (PX  (I)  ,I=1,NP) 

500 

WRITE  (3,  105)  (PY  (I)  ,T=1,NP) 

600 

WRITE( 3,  106)  (PZ  (I), 1=1, HP) 

700 

READ  (1 ,  107)  (LL  (I)  ,I=1,NW) 

300 

WRITE  (3,  108)  (LL  (I)  ,I=1,NW) 

900 

READ(1,1 07)  (LR  (I)  ,1=1,  NR) 

300 

WRITE  (3,  103)  (LP  (I)  ,1=  1 ,  NR) 

100 

READ  ( 1 ,  109)  (RAD  (I)  ,1=1, NR) 

200 

WR ITE  (3 ,  110)  (RAD  (I)  ,1=1, NR) 

300 

101 

FORMAT  (313 , El 4. 7) 

400 

102 

FORM  AT  ( ’ ONW  NP  NR  BK »/3l3 , E 1 4 . 7) 

500 

130 

FORMAT  (10F0.0) 

500 

104 

FORMAT (»  OPX */(1X,8F8.4) ) 

700 

105 

FORMAT  (,0PY,/(1X,8F8.4)  ) 

300 

106 

FORMAT (' 0PZ*/( 1X,8F8.4) ) 

900 

107 

FORMAT  (2013) 

)00 

108 

FORMAT (• OLL */( IX, 1014) ) 

100 

109 

FORMAT  (5E14.7) 

>00 

103 

FORMAT ( *  OLR  * /(IX, 1014) ) 

100 

1 10 

FORMAT  (*0RAD*/(1X,8E0.0)) 

100 

DO  46  1=1,  NR 

500 

RAD2  (I)  =R AD  (I)  *RAD(I) 

500 

46 

CONTINUE 

POO 

J1  =  1 

500 

J2=2 

3900 

1000 

1100 

1200 

4 

1300 

1400 

1500 

1600 

3 

1700 

1800 

1900 

>000 

>100 

>200 

>300 

2400 

>500 

>600 

>700 

>800 

>900 

5000 

2 

noo 

5200 

1300 

1400 

1500 

113 

moo 

noo 

1200 

45 

noo 

noo 

1500 

\  600 

127 

1700 

32 

moo 

1900 

moo 

126 

>100 

1500 

»600 

125 

>700 

1800 

1900 

1000 

100 

.200 

.300 

.400 

.500 

600 

.700 

800 

900 

000 

100 

200 

300 

400 

500 

28 

N1=0 

DO  2  J=1,NP 

IF(LL(J1)-J)  3,4,3 

J2=J 2- 1 

L(J1)=J2 

J1=J1+1 

GO  TO  2 

N  1=N  1+  1 

J3= J- 1 

IF  (  (Nl/2*2—  NT)  .EQ.O)  J2=J2+1 
XI  (N 1)  =.  5*  (PX  ( J)  +  PX  (J3)  ) 

XY  (Nl)  =.5*  (PY  (J)  *PY  (J3)  ) 

XZ  (N 1)  =  .  5*  (PZ  ( J)  +PZ  (J3)  ) 

S1=PX(J)-PX  (J3) 

S2=PY  (J)  -PY  (J3) 

S3=PZ(J)-PZ(J3) 

S4=SQRT  (S1*SUS2*S2  +  S3*S3) 

TX (Nl) =S 1/S4 
TY  (Nl)  =S2/S4 
TZ  (N  1)  =S  3/S4 
AL  <N  1)  =S  4 
CONTINUE 
L(v7l)=J2 

N=J2-2 

CALI  CALZ  (N  ,  N 1 ,  Z) 

NRITE(21, 113)  (Z  (I)  ,1=  1 ,  N) 

FOFMAT  (5E14.7) 

DO  45  1=1, N 
D  (I)  =0. 

CONTINUE 

READ  (1 ,  107)  NSET 
NRITE (3, 127) NSET 
IF  (NSET)  33,33,32 
FORMAT  (»  ONSET ’/I1*) 

DO  26  111=1 , NSET 
U1=(0.  ,1.) 

READ  (1,1 26)  THE,  PHI,  El  (  1)  ,EI(2) 

FORMAT (8E0. 0) 

GUES S= GUESS *R AD  (1)/188.365 
HRITE{3,  125)  THE, PHI, El  (1)  , El  (2) 

FORMAT {fO**********»/’0THETA  PHI 

'  El  (2)  V2F6.0,4E11.4) 

A1=CABS  (El (1) ) **2+CABS (El  (2) ) **2 
A2=BK*BK 
TH=THE*. 0174533 
PH=PHI*. 0174533 
CT=COS  (TH) 

ST=SIN  (TH) 

CP=COS  (PH) 

SP=SIN  (PH) 

S1=CT*CP 

S2=CT*SP 

BK1=BK*ST*CP 

BK2=BK*ST*SP 

BK3=BK*CT 

J1=1 

J2=-2 

DO  27  J= 1 , N 

IF  (L  (J1)  -  J)  29,28,29 

J2=J2*2 


45 


EI(1)  », 


7600 

7700 

7800 

7900 

29 

3000 

3100 

3  200 

30 

3300 

3400 

3500 

3600 

3700 

3800 

31 

3900 

3000 

noo 

32  00 

3  300 

3400 

)500 

)  600 

27 

)000 

)100 

26 

)  200 

33 

3300 

>400 

1  12 

3450 

3500 

3600 

3700 

1200 

1300 

C 

1400 

C 

1500 

C 

1600 

1700 

1800 

1900 

>000 

>100 

>200 

>300 

>400 

>500 

>600 

>700 

>800 

>900 

3000 

4 

3100 

3200 

3 

3300 

3400 

3500 

1600 

3700 

1800 

1900 

3000 

3100 

J1=J1+1  46 

KK=1 

GO  TO  30 
OV{1)=OV  (3) 

OV  (2)  =UV  (4) 

KK=3 

DO  31  M=KK, 4 
J3= J2+M 

XDT=TX  (J3)  *S  1 +TY  (J3)  *S2-TZ(J3)  *ST 

XDP=— TX  ( J3)  *SP  +  TY  (J3)  *CP 

BKB=XX  (J3)  *BK1  +XY  (J3)  *BK2+XZ  (J3)  *BK3 

U7  (M)  =  (X  DT*EI  (  1)  +XDP*EI  (2)  )  *  (COS  (BKR)  +U1*SIN  (BKR)  ) 

CONTINUE 
J3=(J-1)  *4 
J4=J3+1 
J5=J4+ 1 
J6= J5+ 1 
J7=J6+1 

0  ( J)  =T  (  J4)  *UV  (1 )  +T  { J5)  *UV  (2)  +T(J6)  *UV{3)  +T ( J7)  *UV  (4) 

J2=J2+2 

CONTINUE 

WRITE  (21,113)  (U(I),I=1,N) 

CONTINUE 
CONTINUE 
WRITE  (3, 112) 

FORMAT  (1H  , * TY PE  1  TO  STOP,  ELSE  0  AND  RETURN*) 

READ  (1 ,  101)  IFLAG 
IF (IFLAG. EQ.O)  GO  TO  1 
CLOSE  (UNIT=21) 

STOP 

END 

THIS  IS  SUBROUTINE  #1 
SUBROUTINE  CALZ(N,N1,Z) 

COMPLEX  Z  (3600) , PSI (3200) ,H,U1 ,U2,U3,U4,U5,U6 

COMMON  XX  (800)  rXY(800)  ,XZ(800)  ,TX(800)  ,TY(800)  ,TZ(800)  ,AL(800) 
COMMON  T  (16  00)  ,TP(1600)  r  BKr  RAD2  (2  1)  ,  L  ( 1 1)  ,  LR  (2 1) 

DIMENSION  DC  (3200) 

U=  (0. ,  1.) 

PI=3. 141593 
ETA=  376. 7307 
C1=. 125/PI 
C2=. 25/PI 
J1=1 
J2=-  2 

DO  1  J=1,N 

IF  (L  (J1)  -  J)  3,4,3 

J2= J2+2 

J1=J1+1 

J3=  (J-1)  *4 

J4=J3+1 

J5=J4+1 

J6=J5+1 

37= J6* 1 

K4=J2+1 

K5=K4+1 

K6=K5+1 

K7=K6+1 

S1=AL(K4)  +AL  (K5) 


T(J4)=AI  (K4)  *.5*AL(K4)  /SI 

T(J5)=AL  (K5)  *(AL(K4)  +  .  5*AL  (K5)  )  /S  1 

T  (J6)  =AL  (K6)  *  (AL  (K7)  +.  5*AL  (K6)  )  /S2 

T ( J7) =AL (K7) *. 5*AL (K7) /S2 

TP  (J4)  =  AL  (K4)/S1 

TP  ( J5)  =AL  (K5)  /SI 

TP  (J6)=-AL  (K6)  /S2 

TP  (J7)  =-AL(K7)  /S2 

J2= J2+2 

CONTINUE 

IT3=U*BK*  ETA 

U4=- U/BK+ETA 

BK2=BK*BK/2. 

BK3=BK2*BK/3. 

N9=0 
N2=  1 
N0=1 
N3=-  2 

DO  10  NS  =  1 ,  N 

IF (L  (N2) -NS)  12,11,12 

KK=  1 

N3=N3+2 

N2=N2+ 1 

GO  TO  13 

KK=3 

DO  14  NF= 1 , N 1 
N4=N  F+N1 
N5=N4*N1 
N6=N5+N1 
DC (NF) =DC  (N5) 

DC  (N4)  =DC  (N6) 

PSI  (NF)  =PSI  (N5) 

PSI  (N4)=PSI  (N6) 

CONTINUE 
CONTINUE 
DO  15  K=KK, 4 
N7=N3+K 
K1=(K-1)*N1 
N0=  1 

DO  16  NF= 1 , N 1 
IF(NF-LR  (NO)  )  5,6,5 
AA=RAD2 (NO) 

N0=N0+1 

CONTINUE 

N8=NF+K1 

S1=XX(N7)-XX  (NF) 

S2=XY(N7)-XY(NF) 

S3=X7,  (N7)-XZ  (N^) 

R2=S1*SUS2*S2  +  S3*S3  +  AA 
R=SQRT{R2) 

PT=  ABS  (S  1  *rnX  (N7)  *S2*TY(N7)  +S3  +  TZ  (N7)  ) 

RT2=RT*RT 

RH=  (R2-RT2) 

ALP=.5*AI(N7) 

AR= ALP/R 
S  1  =  BK*R 

U2=COS  (S  1)  -U*SIN  (SI) 

IF(AR-.I)  22,22,21 
U2=U2*C1/ALP 


S1=RT-ALP 
52=RT+  ALP 
S3=SQRT (S1*S  1+RH) 

S4=SQRT  (S2+S2  +  RH) 

IF  (Si)  18,18,19 

AI1=AL0G  (  (S2+S4)  *  (-S1  +  S3)  /RH) 

GO  TO  20 

A1 1= ALOG  (  (S2+S4)  /(S1  +  S3)  ) 

AI2=AL  (N7) 

AI3=  (S2*S4-S1*S3  +  RH*AI1)  /2. 

AI4=AI2*  (RH+ALP  +  ALP/3. 4-RT2) 

S3=AI1*R 

S  1= A I 1-BK2*  (AI3-R*  (2.*AI2-S3)  ) 

S2=-BK*(AI2-S3)  *BK3*  ( AI4-3. * AI3*R +R2* ( 3. *AI2-S 3)  ) 

GO  TO  28 

U2=U2*C2/R 

BA=BK*ALP 

BA2=BA*BA 

AR2=AR*AR 

AR3=AR2*AR 

ZR=RT/R 

ZR2=ZR*ZR 

ZR3=ZR2*ZR 

ZR4=ZR3*ZR 

H 1=  (3.-30. *ZR2*35.*ZR4) *AR3/40. 

A1=AR*  (-1.  *3.*ZR2)  /6.  +  (3.-30.*ZR2  +  35.*ZR4)  *AR3/40. 
A0= 1 . +AR*  A  1 

A2=-ZR2/6.-AR2* (1 . -12. *ZR2+ 15. *ZR4) /40. 

A3=AR*  (3.  *ZR2-5.  *ZR4)  /60. 

A4=ZR4/120. 

S1=A0+BA2*(A2+BA2*A4) 

S2=B A*  (A1  +  BA2*A3) 

PSI  (N8) =  U2*  (S1+U+S2) 

DC  (N8)  =TX  (NF)  *TX  (N7)  +TY  (NF)  *TY  (N7)  +TZ  (NF)  *TZ  (N7) 

CONTINUE 

CONTINUE 

N3=N3+2 

J3=  (NS-1) *4 

J7=-  2 

J9=  1 

DO  25  NF= 1 , N 

Jl  = (NF- 1 ) *4 

IF  (L  (J9)  -NF)  26,27,26 

J9=J9+ 1 

J7= J7+2 

N9=N9+ 1 

05=0. 

U6=0. 

J5=0 

DO  23  JS= 1 , 4 
J4=J3+JS 
J8=J5+J7 
DO  24  JF= 1,4 
J6=J8+ JF 
J2=JUJF 

U5=T  (J2)  *T  ( J4)  *DC  (J6)  *PSI  ( J6)  +U5 

06=TP(J2)  *TP(J4)  *PSI  (J6)  +06 

CONTINUE 

J5=J5+N1 

CONTINUE 
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10 

Z(N9)=U5*U3+U6*U4 

J7=J7+2 

IF  (N9.  GE-  N)  RETURN 

CONTINUE 

CONTINUE 

RETURN 

100 

110 

$ 

$ 

END 

SUBROUTINE  HELIX <X, I, Z , NP) 

REAL  X  (1)  ,Y  (1)  ,Z  (1)  ,  RAD,  PITCH,  PHI,  DPHI 
INTEGER  TURNS 

READ  (1 , 100)  RAD, PITCH, TURNS 

WRITE (3, 110) RAD, PITCH, TURNS 

FORMAT  (2F0.  0,10) 

FORM  AT (1 H  , *RADIUS=*,E14.7/ 

1 H  , *  PITCH  =* , E14. 7/ 

1 H  ,  '  NO  OF  TURNS=’  ,17) 
PITCH=PITCH/6. 2831853 

PHI=0. 0 

DPHI=6. 2831853*TURNS/FLOAT (NP-1) 

DO  10  1=1, NP 

X  (I) =RAD*COS (PHI) 

Y (I) =RAD*SIN{PHI) 

Z  (I)  =PITCH*PHI 

PHI=PHI+DPHI 

10 

CONTINUE 

RETURN 

END 
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XI.  MAIN  PROGRAM  SEGMENT  AND  SUBROUTINES  FOR  THE  PLANAR 

ARRAY  PROBLEMS  (  SINGLE  EXPANSION  FUNCTION  PER  ANTENNA 
ELEMENT  ) 


Subroutine  SHAPE  zerorize  the  corners  i.e.  it  zerorize 
the  phantom  antenna  element  currents.  Subroutine  PATTRN 
computes  the  array  factor  for  the  given  set  of  array 
currents.  The  main  program  segment  solves  the  planar  array 
problems  using  the  solution  routine  SOLVE  and  other  routines 
not  listed  in  this  section.  They  are  listed  in  section  III 
since  these  routines  are  common  to  both  single  and  multiple 
expansion  solution  programs. 

The  flow  chart  for  the  main  program  segment  is  as 
follows . 


51 


ip  be  solved?. 


STOP 
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PROGR  Art  TO  COMPUTE  MUTUAL,  IMAGE  AND  TOTAL  IMPEDANCES 
LOGICAL  FXCITE 

INTEGER  CODE, FLAG, ECSPTE,FLAG2 

COMPLEX  MUTUAL, IMAGE, ZMNG,V (1849)  ,Z  <16384) ,CZERO, Y { 1849) 
READ  IN  N ROWS  -  NUMBER  OF  ROWS  OF  ANTENNA  ELEMENTS 
NCOLS  -  NUMBER  OF  COLUMNS 
DX  -  SEEEP ATION  EETNEEN  ELEMENTS  IN  X  DIRECTION 
DY  -  DISTANCE  OF  ELEMENTS  FROM  GROUND  PLANE 
DZ  -  SEFERATION  EETWEEN  ELEMENTS  IN  Z  DIRECTION 
FXCITE=-TRUE. 

TWOPI=6. 2831853 
CZERO=  (0-0, 0.0) 

OEEN  (UNIT=21, FILE='P ATTRN.DAT') 

READ(1,  100)  NSIDE, CODE, D X , D Y , DZ , WLNGTH , HAD 
N T=NSIDE*  3-2 
NPO  WS  =  NT 
NCOLS=NT 

CALL  SIZE  (NCOLS, NROWS,LM, MM, L,M,N) 

DC  1  1=1, N 
Z  (I) =CZERO 
CONTINUE 
DX=DX*TWOPI 
D Y=DY*TWO  PI 
DZ=DZ*TKOPI 
WLNGTH=WLNGTH*TWOEI 
RAD=R AD*TWOPI 
NILEM=NBOWS*NCCLS 
ZFLEM=C. 0 
D Y2SQ=4. 0*DY*EY 
R AD2=EAD*RAD 
WLBY2=WINGTH/2- 0 
POSPTR= 1 

DO  20  1=1,2*  NT- 1 
XCUR=  (I- 1 ) *DX 
ZCUF=  (6*NSIDE— 5-I)*DZ 
DO  10  J=  1 , 2*  NT—  1 
RR=XCUR*XCUR+EY2SQ 

IMAGE=ZMNG  (ZELEM,ZCUE, WLBY2, BLEY2 ,RR) 

RE=RR-DY2SQ 

IF  (RR-EQ-0.0)  RR=RB  +  RAD2 
MUTUAL=ZMNG  (Z EL EM ,ZCUR, WLEY2 , WLB Y 2 , RR) 

Z  (POSPTR) = MUTUAL— IMAGE 
POSET F=POS  PTR+ 1 
XCUR=  XCUR—  DX 
ZCUR=ZCUR— DZ 
CONTINUE 

POSPTB=POSPTR+L— NT— NT+ 1 
CONTINUE 

R E  AD  (  1 , 100)  CODE, El AG 

IF (CODE. EQ. 1 )  CALI  VCLTU  ( NCOLS , NRCWS, V) 

IF  (CODE.  EQ.  2)  CALL  VCLTC  (NCOLS, NROWS, V) 

IF (CODE. EQ. 3)  CALI  TA P  (  NCGLS , NBO WS , DX , DZ , V) 

IF  (CODE.  EQ.  4)  CALI  VCLTK  (NCOLS, NROWS,  V) 

IF (CODE. EQ. 5)  CALL  STEER  ( NCOLS , NBO WS , DX , DZ , V) 

IF  (CODE.  EQ. 6)  CALL  PHASE  (NCOLS, NROWS,  V) 

CALL  SOLVE  (Z , V, Y , NCOLS , NROWS , NELEM ,LM , M H, L, M, N, FXCITE) 
FXCITE=. FALSE. 

WRITE  (3,300) 

REACH, 100)  FIAG2 
IF (FLAG2. NE. 0)  GO  TO  25 
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CALI  E ATT BN  ( NCOLS , N RO WS, N EL  EM , Y,DX,DZ) 

CALL  SHAPE  (V,NSIDF, NT) 

CALL  PATTRN  ( NCOLS , NRO KS , NELEM , V , DX , DZ) 

CCNTINUE 

IF  (FLAG. EQ. 0)  GC  TC  30 

STOP 

100 

200 

300 

FORMAT (210, 5F0.  0,10) 

FC  RM AT  ( 1 H  ,  4  E  1 2.  4 ) 

FORMATC1H  ,'TYPE  0  TO  PRINT  ARRAY  FACTOR,  ELSE  1») 
END 

c 

c 

RGUTINE  TO  COMPUTE  THE  FAR  FIELD  PATTERN 

SUBROUTINE  PATTRN  (10 , MO , NO, Y  , D X , DZ) 

INTEGER  LC , MO , NC , FTB 

REAL  DX,DZ,ALPHA,CCSTHE, MAGE, DANG, T1,T2,T3, DEL 
COMPLEX  Y  (1)  , AF,XPFASE,ZPHASE,EXPH,DZPH,CNPLX 

REAL  PAT  (361) 

WRITE  (3,1000) 

READ{1,  1100)  NPTS  ,  D  A  NG 

AIPHA=0. 0 

DO  300  1.1=1,  NETS 

COSTHE=COS (ALPHA*. 1745329 3E-0 1) 

I 3=DX*COSTHE 

T  1  =  COS  (T3) 

T  2=SI  N  (T  3) 

DXPH=CMPIX  (T  1  ,T2) 

DZPH=DXPH 

A  F=  (0.0, 0.0) 

ZPHASE= (1.0, 0.0) 

DC  200  1=1, MO 

PTR=  (I-  1)  *LO 

XPH AS  E=  (1. 0,0.0) 

DO  100  .1=1,  LO 

AE=AF  +  Y  (PTR  +  J) *XPH AS E*Z PH AS E 

XPHASF=XPHASE*DXPH 

100 

CCNTINUE 

ZPHASE=ZPHASE*DZPH 

200 

CONTINUE 

M AGE=CABS (AF) 

PAT  (II)  =MAGE 

WPITE (3, 1200)  ALPHA, MAGE, AF 

ALPHA =ALPHA+EANG 

300 

CONTINUE 

1000 

1  100 
1200 
1300 

SPITE  (21, 1300)  (EAT  (II)  ,11=1,  NPTS) 

FORM  AT  ( 1 H  , 'NPTS  AND  ANGLE  INCREMENT?') 

FORMAT  (10, EO.  0) 

FORMAT (  1H  ,F8.1,3F12.4) 

FORMAT  (8E11.4) 

RETURN 

END 

c 

c 

ROUTINE  TO  ZEROIZE  CORNERS 

10 

SUBROUTINE  S HAPE  { Y , NS , LO) 

COMPLEX  Y  (1)  , CZ  E  RC 

CZERO=  (0.0, 0.0) 

DO  30  1=1, NS— 1 

IPTR=  (I-  1)  *LO 

DO  10  J= 1 , NS-I 

Y (IPTR+J) =CZERO 

CCNTINUE 
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12100 

DO  20  J=2*NS+I-1#3*NS-2 

12200 

Y  (IPTB  + J) =CZ  ERO 

12300 

20 

CONTINUE 

12400 

30 

CONTINUE 

12500 

DO  60  I=2*NS,3*NS-2 

12600 

IPTR=  (1-1)  *10 

12700 

DO  40  J=1,I-2*NS+1 

12800 

Y  (IPTR+J) =CZERC 

12900 

40 

CONTINUE 

13000 

DO  50  J=5*NS-I-2,  3*NS-2 

13100 

Y  (IPTE  +  J)=CZEBO 

13200 

50 

CONTINUE 

13300 

60 

CONTINUE 

13400 

RETURN 

13500 

END 

00  050 

SUBROUTINE  SOLVE (A, B, LO , MO, NO, LM , MM ,L , H, N, EXCITE)  °5 

00100 

C 

BOUTINE  TO  SOLVE  THE  MATRIX  EQUATION  A  X  =  B 

00150 

LOGICAL  EXCITE 

00200 

COMPLEX  CTEMP,CZERO,  A  ( 1 )  ,X(4096)  ,V  (4096)  ,B(1)  ,Y(7  396) 

00300 

INTEGER  FLAG 1 , FLAG 2, COUNT 

00  400 

REAL  CHNAVG,CHNMAX, CHANGE 

00600 

CZERO={0.  0,0.0) 

01 100 

c 

ZEROIZE  Y  AND  V  (EXPENDED  B) 

01300 

DO  10  1=1, N 

01400 

V  (I)  =CZERO 

01500 

10 

CONTINUE 

01600 

DO  20  1=1, NO 

01700 

Y  (I)  =CZERO 

01800 

20 

CONTINUE 

01  900 

IF  (.  NOT.  FXCITE)  GO  TO  65 

02000 

DO  40  1=1, MO 

02100 

IPTR=  (1-1)  *L 

02300 

JPTR=IPTR+LO+LO-1 

02400 

DO  30  J=IPTR+ 1, IPTR+LO- 1 

02500 

A  (JPTR)  =A  (J) 

02600 

JPTR= JPTR- 1 

02700 

30 

CONTINUE 

02800 

40 

CONTINUE 

03000 

c 

FILL  UP  A  ARRAY  AND  V  ARRAY 

03100 

DO  60  1=1, MO-1 

03200 

IPTR=  (1-1)  *L 

03300 

JPTR=  (MO+MO-I-1)*! 

03400 

DO  50  J= 1 , LO+LO- 1 

03500 

A  (JPTR+J)  =  A  (IPTR+J) 

03600 

50 

CONTINUE 

03  700 

60 

CONTINUE 

03900 

CALL  TWODF  (A,N,M,L,MM,LM) 

04000 

65 

CONTINUE 

04  005 

CALL  ITFR  (A,B,X,V,Y, LO, MO,NO,LM,MM,L,M,N) 

04010 

RETURN 

04015 

END 

04  02  0 

SUBROUTINE  ITER  (A,  B,  X,  V ,  Y,LO ,  MO ,  NO  ,LN  , M M,L,  M  ,  N) 

04  025 

INTEGER  COUNT, FLAG1,FLAG2 

04  03  0 

COMPLEX  C7ERQ,CTEMP,  A  {  1)  ,B(1)  ,X(1)  ,V(1)  ,Y(1) 

04035 

REAL  CHNAVG,CHNMAX, CHANGE 

04040 

CCUNT=0 

04045 

CZERO=(0.  0,0.  0) 

04100 

70 

JPTR= 1 

04200 

COUNT=COUNT+ 1 

04300 

DO  90  1=1, MO 

04400 

IPTR  =  L* (MO-2+1)  +LO 

04500 

DO  80  J=1 , LO 

04600 

V (IPTR) =B  (JPTR) 

04700 

JPTR= JPTR+  1 

04800 

IPTB=IPTR+  1 

04  900 

80 

CONTINUE 

05000 

90 

CONTINUE 

05200 

c 

FIND  V  TRANSFORMED  AND  COMPUTE  X  TRANSFORMED 

05300 

CALL  TWO DP (V,N,M,L,MM,LM) 

05400 

DO  100  1=1,  N 

05500 

X  (I)  =V  (I)  /A  (I) 

05600 

100 

CONTINUE 

05700 

c 

GET  X  FROM  X  TRANSFORMED 

05800 

CALL  IT  WODF  (X,N,M,L»MM,LM) 

06100 

c 

TRUNCATE  X  AND  SAVE  X  AFTER  COMPUTING  THE  CONVERGENCE  CRITERION 

06200 

06300 

064  00 

06500 

06600 

06700 

06  800 

06900 

07  000 

07100 

07200 

07300 

07400 

07500 

07600 

07700 

07800 

07  900 

08000 

08100 

08  200 

08300 

08400 

08500 

08600 

08700 

08  800 

08900 

09  000 

09100 

09200 

09  300 

09400 

09500 

09600 

09700 

09  8  0  0 

09900 

10000 

10100 

10200 

10300 

10400 

10500 

10600 

10700 

10800 

10900 

11000 

11 100 

11200 

11300 

11400 

11500 

11600 

11700 

11800 

11900 

12000 

12100 
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CHNAVG=0. 0 
CHNMAX=0. 0 
DO  120  1=1,  N 
IPTH=I/L 
JPTB=I— IPTR*L 

IF(JPTR.LE.LO.AND. JPTR.NE.O. AND. IPTR. LT. HO)  GO  TO  110 
X  (I)  =CZERO 
GO  TO  120 

110  IPTR=IPTR*LO+ JPTR 
CTEMP=X  (I) 

CHANGE=CABS  (CTEMP-Y  (IPTR)  )  /CABS  (CTEMP) 

Y (IPTR) =CTEMP 
CHNAVG=CHNAVG+CHANGE 
IF  (CHANGE. GT.CHNM  AX)  CHNMAX=CHANGE 
120  CONTINOE 

CHNAV G= (CHNA VG* 100.0) /FLOAT  (NO) 

CHNMAX=CHNMAX* 1 00.0 
WRITE  (3,1200)  CHNAVG,CHNMAX,COONT 
C  FIND  THE  TRANSFORM  OF  TRONCATED  X 

CALL  TWODF  (X ,  N,  M ,  L,  M M,  LM) 

C  COMPOTE  V  TRANSFORMED 

DO  130  1=1, N 
V(I)=A(I)  *X(I) 

130  CONTINUE 

C  GET  V  FROM  V  TRANSFORMED 

CALL  ITWODF (V,N,M,L,HM,LM) 

C  COMPOTE  THE  ERROR  CRITERION 

CHNAVG=0. 0 
CHNMAX=0. 0 
JPTR=  1 

DO  150  1=1, MO 
IPTR=L* (MO-2+I) +LO 
DO  140  J= 1 , LO 
CTEMP=B  (JPTR) 

CHANGE=CABS (CTEMP-V (IPTR) ) /CABS (CTEMP) 

C  HNA VG=CHNAVG+CHANGE 
IF (CHANGE. GT.CHNM AX)  CHNH AX=CHANGE 
IPTR=IPTR+  1 
JPTR= JPTR+ 1 
140  CONTINUE 
150  CONTINOE 

C  ASK  WHETHER  OR  NOT  TO  STOP  AFTER  REPORTING  *  FIELD  ERROR 

CHNAVG=  (CHNAVG*  100. 0) /FLOAT (NO) 

CHNMAX=CHNMAX* 100.0 
WRITE  (3, 1300)CHNMAX, CHNAVG 
READ(1,  1100)FLAG1 
IF  (FLAG  1. EQ. 0)  GO  TO  70 
WRITE(3,  1400)  (T  (I)  ,1=1, NO) 

C  ASK  IF  FIELD  SHOULD  BE  PRINTED  OOT  ALSO 

WRITE  (3,  1500) 

READ (1,1 100) FLAG2 
IF (FLAG2. HE. 0) RETURN 
WRITE  (3, 1600)  (V  (I)  ,  1=  1 ,  N) 

RETURN 

1000  FORMAT (10E0.  0) 

1100  FORMAT  (410) 

1200  FORMAT  (1H  , ' AVG  CURRENT  CHANGE= ' , El 4. 7, •  * •/ 

$  1H  , *  MA X  CURRENT  CHANGE=f , E 1 4. 7 , *  %•/ 

$  1 H  , ’AFTER', 14,'  ITERATIONS'/) 

1300  FORMAT  (1H  ,»MAX  FIELD  ERROR  =  '#E15.7,'  %*/ 


12200 

12300 

12400 

12500 

12600 

12700 
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$  1 H  ,  »AVG  FIELD  ERROR  =  ',E15.7,'  *'/ 

$  1H  , ’CONTINUE  ITERATIONS?  0  FOR  YES,  1  FOR  NO,  AND  RETURN'/) 

1400  FORMAT  ( 1H  ,* CURRENTS '//{ 1H  ,10E11.4)) 

1500  FORMAT {1H  PRINT  FIELDS?  0  FOP  YES,  1  FOR  NO,  THEN  RETURN'/) 
1600  FORMAT  (1H  , 'RESULTANT  FTELDS'//(1H  ,10E11.4)) 

END 
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III.  MAIN  PROGRAM  SEGMENT  AND  SUBROUTINES  FOR  THE  PLANAR 

ARRAY  PROBLEMS  (  THREE  EXPANSION  FUNCTIONS  PER  ANTENNA 

ELEMENT  ) 

The  difference  between  the  main  program  segment  for 
three  expansion  functions  per  antenna  element  and  the  main 
program  segment  for  one  expansion  function  per  antenna 
element  is  that  there  are  five  distinct  mutual  impedance 
vectors  z  to  be  computed  for  the  three  expansion  functions 
per  antenna  element  solution.  Therefore  the  computation  of 
each  Z  is  done  by  the  seperate  routine  FILLZ.  Similarly, 
MSOLVE  and  SOLVE  routines  differ  mainly  in  that  there  are 
three  distinct  vectors  each  for  the  generalized  voltage  V 
and  J  to  be  computed  by  MSOLVE.  SOLVE  on  the  other  hand, 
computes  only  one  vector  each  of  V  and  J. 

Routine  ZMNG  computes  the  mutual  impedance  between  two 
parallel  segments  of  thin  wire  dipoles  of  same  length  which 
may  or  may  not  be  offset  from  each  other  along  one  or  more 
axes.  SICI,  VOLTU ,  VOLTC  TAP,  VOLTK ,  and  PHASE  are  all  from 
[9]  and  hence  no  description  of  them  will  be  given  here. 
FFT,  IFFT,  TWODF,  and  ITWODF  are  fast  fourier  transform  and 
inverse  transform  routines  for  one  and  two  dimensional 
discrete  fourier  transforms,  respectively.  PATT3E  is  the 
routine  to  compute  the  array  factor  from  the  current 
distribution  solutions  obtained  by  using  three  expansion 
functions  per  antenna  element. 
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The  last  main  program  segment  computes  the  array 
factors  from  the  current  distribution  solutions  obtained  by 
using  three  expansion  functions  per  antenna  element  and  from 
the  ideal  solutions  which  ignore  the  mutual  coupling  between 
antenna  elements. 


PROGRAM  FOR  TRIANGULAR  PATTERN  WITH  THREE  EXPANSIONS  PER 
ELEMENT 

LOGICAL  FXCITE 

INTEGER  CODE, EL AG, NRCWS , NCGLS, NE1EM 
REAL  TWCPI 

COMPLEX  V 2  (4  84)  ,ZA  (4  096)  , ZB  (4096)  ,ZC{4096)  ,ZD  (4096)  ,ZE(4096) 
COMPLEX  XI  (4096)  ,X2  (4096)  ,X3(4096) 

FXCITE-=.  TRUE. 

TKOPI=6. 2831853 
CZEBO=  (0.  0,0.  0) 

OEEN (UNIT=21 ,FILE=' MXCUR- CAT') 

R  EAD (  1,  100)  N SIDE, DX,DY,DZ,WLNGTH, RAD 

NHOWS=NSIBE*  3-2 

NCOLS=NROHS 

CALL  SIZE  (NCOLS, NRCWS,LM,MM,L,H,N) 

DX=DX*TWOPI 

DY=DY*TWOPI 

DZ=DZ*TWOPI 

WI.BY2=WLNGTH*TWO  PI/4.0 

ZFLEM=WLFY2 

RAD=RAD*TWOPI 

N EL EM=NROWS* NCOLS 

DY2SQ=4. 0*  DY*EY 

R  AD2=BA  D*R  AD 

ZCF  F=  0.  0 

CALL  FILLZ (ZA , ZOFF , DX , DZ , DY2SQ, WLBY2 , R AD2, NROWS ,L, N) 
ZOFF=-ZELEM 

CALL  FILLZ (ZB , ZCFF , DX , DZ , DY2SQ, WLBY2, B A D2, NRC8S, L, N) 

ZOF  F=—  ZELEM-ZELEM 

CALL  FILLZ (ZC,ZOFF, DX , DZ ,DY2 SC , WLBY2, EAD2, NRCWS , L, N) 

CALL  EIILZ  (ZD,Z  EL  EM,  DX,  DZ  ,  DY2SQ,  WLBY2  ,R  AD  2,  NROWS  ,L  ,N) 
ZOFF=ZELEM+ZELEM 

CALI  FILLZ (ZE,ZCFF, EX, DZ,DY2SQ,WLBY2,RAD2, NRCWS,  L,N) 

R  EAD  {  1,  200)  CODE, FLAG 

IF  (CODE.  EQ.  1 )  CALL  VCLTU  (  NCOLS  ,NRO  HS,  V 2) 

IF  (CODE. EQ. 2)  CALL  VCLTC  (NCOLS, N ROWS, V2) 

IF  (CODE. EQ. 3 )  CALI  T AP ( NCOLS , NROWS , DX , DZ , V2) 

IF (CODE. EQ. 4)  CALL  VCLTK (NCOLS , NRCWS, V2) 

IF  (CODE. EQ. 5)  CALL  STEER (NCOLS , NROWS , DX , DZ , V2) 

IF (CODE. EQ. 6)  CALL  PH ASE  { NCOLS , NROWS , V 2) 

CALL  MSOLVE (Z A, Z B, ZC,ZD, Z E , V2 , X 1 , X 2 ,X 3, NCOLS , NRO WS , 
NELEM,LM,MM,L,M,N, FXCITE) 

FXCITE=. FALSE. 

IF  (FLAG- EQ. 0)  GO  TC  10 
STOP 

FORMAT (I0,5F0.0,IC) 

FORMAT  (210) 

END 

ROUTINE  TO  FILL  THE  Z  MATRIX  FCR  TRIANGULAR  ARRAYS 
SUBROUTINE  FILLZ  (Z , ZCFF , DX , DZ , EY2SQ , WLBY2 ,R AD2 , NRO WS,L, N) 
INTEGER  POSPTR, NRCWS, L, N 

REAL  DX,CZ ,DY2SQ„ZELEM,XCOR,ZCUR,RR,WLBY2 
COMPLEX  Z(1)  ,  MUTUAL,  IM  AGE  ,  CZEBO  ,  Z  M  NG 
CZ  ERO=  (0.0,0.  0) 

ZELEM=0. 0 
PCSPTR= 1 
DO  10  1=1, N 
Z  (I) =CZERC 
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10  CONTINUE 

NTOTAl=NROWS*NRCW£-1 

DO  30  1=1, NTOTAL 

XCUR  =  (1-1)  *DX 

ZCUR=  (NTOTAL-I) +DZ+ZCFF 

DO  20  J  =  1 , NTOTAI 

RR=XCUE*XCUR+DY2SQ 

IMAGE=ZMNG  (ZEIEM , ZC UR, WLE Y2 , EB) 

BR=  RR-DY2SQ 

IF  (RR-EQ. 0.0)  RR=RE+RAD2 
MUTTJAL=ZMNG  (Z  E L  E  M ,  ZC U R ,  WL  E Y 2  ,  R  R) 

Z  (POSPTR) =  MUT  UA  L-I M  AGE 
POSPTR=FCSPTR+ 1 
XCUR  =  XCUR— DX 
ZCUR=ZCUF-DZ 
20  CONTINUE 

PO SPT 5= POSPTR+L- NTOTAL 
30  CONTINUE 
RETURN 
END 
C 

C  ROUTINE  TO  SOLVE  TEE  CONVOLUTICN  EQUATIONS 

SUBROUTINE  HS01VE  (ZA,ZB,ZC,ZE,ZE,B,Xl,X2,X3,LO,HO,NO, 

$  LM,NM,L,M,K, FXCITE) 

INTEGER  COUNT 
LOGICAL  FXCITE 

CCMPLEX  CTEHP1,CTEMP2,CTEMP3,ZA  (1)  ,ZB  (  1)  ,ZC  (1)  ,ZD  ( 1)  ,ZE  (1)  , 

$  B  ( 1)  ,X1  (1)  ,X2  (1)  ,X3  (1)  ,V1  (4096)  ,V2(4096)  ,73(4096)  , 

$  T1 ,T2,T3,T4,T5,T6,T7,T8,T9,T10,Y1 (484) ,Y2 (484) ,Y3 (484) , 

$  CZERO 

CZERC=  (0. 0,0. 0) 

DO  10  1=1, NO 
Y 1  (I) =CZERC 
Y 2  (I) =CZERO 
Y3  (I) =CZERO 
10  CONTINUE 

NS=  (LC  +  2)/3 

EFACT=3 . 0  + FLO  AT  (NO-2*(NS-1)*NS) 

DO  20  1=1,  N 

V  1  (I) =CZERO 
V2  (I) =CZERO 
V 3  (I) =CZERO 

20  CONTINUE 

IF  (.NOT. FXCITE)  GC  TO  30 
CALL  TWOEF  (Z  A ,  N  ,  M,  L ,  M M  ,LM) 

CALL  TWODF  (ZB,N,M,I,Mt5,LN) 

CALL  TWODF  (ZC,N,fi,L,  MM,LM) 

CALL  TWODF  (ZD,N,f1,L,N?1,LM) 

CALL  TWODF  (ZE,N,  M,I,  F«,L!1) 

30  ICENTR=L*  (EO- 1) +LC- 1 
CGUNT=  1 
40  CONTINUE 

DO  60  1=1 , NS- 1 
IPTR=  (1-1)  *L+ICENTR 
JPTE=  (1-1)  *LO 
DO  50  J=NS— I  +  1,2*NS  +  I—  2 

V  1  (IPTR  +  J) =CZ  ERC 

V2  (IPTR  +  J) =B (JPTR+J) 

V3  (IPTR  +  J)=CZERG 
50  CONTINUE 
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11600 

60 

CONTINUE 

11700 

DO  80  I=NS , 2*  NS- 1 

11  800 

IPTR=  <I-1)*L  +  ICENTE 

11900 

JPTR  =  <1-  1)  *LO 

12000 

DO  70  J=1,3*NS-2 

12100 

V  1  (IPTR+ J) =CZEBO 

12200 

V 2 (IPTR+J) =  B  (  JPTR  +J) 

12300 

V3 (IPTR+J) =CZERO 

12400 

70 

CONTINUE 

12500 

80 

CONTINUE 

12600 

DC  100  I=2*NS,3*NS-2 

12700 

IPTR=  (1-1)  +L+ICENTR 

12800 

JPTR=  (1-1) *10 

12900 

DO  90  J=I+2-2*NS,5*NS— 1-3 

13000 

VI  (IPTR  +  J) =CZERC 

13100 

V 2  (IPTR  +J ) =B  ( JPTR+J) 

13200 

V  3 (IPTR+J) =CZERC 

13300 

90 

CONTINUE 

13400 

100 

CONTINUE 

13450 

105 

CONTINUE 

13500 

C 

FIND  V 1  ,  V  2  ,V  3  TRANSFORMED  AND  COMPUTE  Xl,X2,X3 

13600 

CALL  TWODF  (V  1 ,  N  ,  fl  ,1 ,  M  M  ,  LM) 

13700 

CALL  TWODF  (V 2  ,  N,  M , I,  MM  ,  LM) 

13800 

CALL  TWODF(V3,N,M,I,MM,LM) 

13900 

DC  110  1=1,  N 

14000 

T 1  =  ZD (I) /ZA  (I) 

14100 

T2=ZE  (I)/ZA  (I) 

14200 

T3=ZA  (I)  -T  1*ZB  (I) 

14300 

T4=ZB  (I)  -T  1*ZC  (I) 

14400 

T 5=ZA (I)  — T2*ZC  (I) 

14500 

T6=ZD  (I)  -T2*ZB  (I) 

14600 

T7=V2  (I)-V1  (I)  *T  1 

14700 

T8=V3  (I)-V1  (I)  *T2 

14800 

T9=T5-T4*T6/T3 

14900 

T10=T8-T7*T6/T3 

15000 

X  3 (I) =T 10/T9 

15100 

X2  (I)=  (T7-T4  +  X3  (I)  )/T3 

15200 

XI  (I)  =  (VI  (I)- ZB  (I)  *X2  (I)-ZC(I)  *X3  (I))/ZA  (I) 

15300 

110 

CONTINUE 

15400 

CALL  ITWODF (XI, N, M,L,MM,LM) 

15500 

CALL  ITWODF  (X2,N,M,L,  MM, LM) 

15600 

CALL  ITWODF  (X3,N,M,L, MM, LM) 

15700 

DC  140  1= 1 , NS-  1 

15800 

IPTR=  (1-1)  *1 

15900 

DC  120  J= 1 , NS— I 

16000 

X  KIPTR+J)  =CZFRO 

16100 

X2  (IPTR+J) =CZFGC 

16200 

X3  (IPTR  +  J) =CZERO 

16300 

120 

CONTINUE 

16400 

DO  130  J=2*NS+I-1 , 3*NS-2 

16500 

X  1  (IPTR+J) =CZERC 

16600 

X2  (IPTR+J) =CZERC 

16700 

X3  (IPTR+J) =CZERO 

16800 

130 

CONTINUE 

16900 

140 

CONTINUE 

17000 

DO  170  I=2*NS , 3*NS-2 

17100 

IPTR=  (1-1)  *L 

17200 

DO  150  J=  1 ,1-2*  NS+ 1 

17300 

XI  (IPTR+J) =CZERC 

17400 

X2  (IPTR  +  J) =CZERC 
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17500 

X  3  (IPTR+J)=CZERC 

17600 

150 

CONTINUE 

17700 

DO  160  J=5*NS— 1-2, 3* NS- 2 

17800 

X  1  (IPTR+J) =C  Z  E  R  C 

17900 

X2  {IPTR  +  J ) =CZ  EEC 

18000 

X3  (IPTR+J) =CZERC 

18100 

160 

CONTINUE 

18200 

170 

CONTINUE 

18300 

C 

18400 

C 

COMPUTE  THE  CRITERIC NS  AND  REPCBT 

18500 

C  HN  A  V  G=  0. 0 

18600 

CHN  M  A  X=  0.  0 

18700 

DC  190  1=1,  N 

18800 

IPTR= I/L 

18900 

JPTR=I-IPTR*I 

19000 

IF (JPTR.LE.LO.AND. JPTR. NE.O. ANE.IPTR. IT. MC)  GOTO  180 

19100 

X 1  (I) =CZERO 

19200 

X2  (I)  =CZERO 

19300 

X 3  (I)  =CZE RO 

19400 

GO  TO  190 

19500 

180 

CONTINUE 

19600 

IPTR=IPTR*LO+JFTR 

19700 

CT  EM  P 1  =  X 1  (I) 

19800 

CTEMP2=X2  (I) 

19900 

CTEMP3=X3  {I) 

20000 

IF  (CTEMP1 - EQ.CZERO)  GO  TC  190 

20100 

CHNGE1=CABS  (CTEMPl-YI  (IPTR)  )  /CAES  (CTEMP1) 

20200 

CHNGE2=CAB5(CTEMP2-Y2  (IPTR)  )/CABS  (CTEEP2) 

20300 

CHNGE3=CABS  (CTEMP 3-Y  3  ( IPTR)  )  /CAES  (CTEHP  3) 

20400 

Y 1 (IPTR) =CTEMP 1 

20500 

Y2 (IPTR) =CTEMP2 

20600 

Y 3  (IPTR)  =CTEPIP3 

20700 

CHNAVG=CHNAVG+CHNGE  1 +CHNG E2+CHNGE3 

20800 

IF  (CHNGE1-GT.CHNMAX)  CHNMAX=CHNGE  1 

20900 

IF  (CHNGE2.GT.CHNMAX)  CHNMAX=CHNGE2 

21000 

IF (CHNGE3.GT.CHNMAX)  CHNMAX=CHNGE3 

21  100 

190 

CONTINUE 

21200 

CENAVG=  (CHNAVG*  10  0. 0)  /EFACT 

21300 

C H N M A X=  (CHN MAX*  100.  0) 

21  400 

WRITE  (3,  1200)  CHNAVG,CHNMAX,CCUNT 

21500 

C 

21600 

C 

GET  TEE  TRANSFORM  OF  THE  CURRENTS 

21700 

CALI  TWODF  (X1,N,M,I,MM,LM) 

21800 

CALL  TWODF (X2,N,M,L,KM,LM) 

21900 

CALL  TWODF (X3,N,M,L,MM,LM) 

22000 

DO  200  1=1, N 

22100 

VI  (I)  =ZA  (I)*X1  (I)  +ZB  (I)  *X2(I)  +ZC(I)  *X3  (I) 

22200 

V 2  ( I)  =ZD  (I)  * X  1  (I)  +ZA  (I)  *X2  (I)  +ZB  (I)  *X3  (I) 

22300 

V3  (I)  =  ZE  (I)  *X1  (I)  +ZD  (I)  *X2  (I)  +ZA  (I)  *X3  (I) 

22400 

200 

CONTINUE 

22500 

c 

22600 

c 

GET  TEE  FIELD  AND  COMPUTE  %  ERRORS 

22700 

CALL  ITVODF(Vl,N,fl,L,nfl,Lfl) 

22800 

CALL  ITWODF (V2,N,M,1,MM,LM) 

22900 

CALL  ITKODF{V3,N,8,L,MM,LM) 

23000 

CHNAVG=0. 0 

23100 

CHNMAX=0.  0 

23200 

DO  220  I  =  1 , N S—  1 

23300 

I PTR=  (1-1)  *L  +  ICENTE 

23400 

JPTR=  (1-1)  *LO 

DO  210  J=NS-I+1,2*NS+I-2 
CTEMP1=B  (JPTR+J) 

CHANGE=CABS  (CTEMP1-V2  (IPTR  +  J)  )  /CABS  (CTEMP1) 

IF  (CHANGE. GT.CHNMAX)  CHBH AX=CHANGE 
C  HN  A  VG=CHN A VG+CHA  NGE 
V2  (IPTR+J)=CTEME1 

V  1  (IPTR+J) =CZEBC 
V3  (IPTR+J) =CZ  EEC 

210  CONTINUE 
220  CONTINUE 

DO  240  I=NS , NS*2—  1 
IPTR=  (1-1)  *L+ICENTE 
J PTR=  (1-1)  *10 
DO  230  J=1,3*NS-2 
CTEHP 1=B ( JPTR+J) 

CHANGE=CABS  (CTEM P  1- V 2 ( IPTR + J)  )  /CABS  (CTEMP1) 

IF  (CHANGE.  GT.CHNflAX)  CHNMA  X=CHANGE 
CHNAVG=CHN A VG+CHA NGE 
V 1 (IPTR+J) =CZEBG 
V2  (IPTR+J)=CTEMP1 

V  3 (IPTR+J) =CZERC 
230  CONTINUE 

240  CONTINUE 

DO  260  1=2* NS, 3* NS— 2 
IPTR=  (1-1)  *L  +  ICENTB 
J  FT E=  (1-1)  *LO 
DO  250  J=I+2-2*NS,5*NS-I-3 
CTEMP1=E  (JPTR+J) 

CHAN  GE=CAB  S  (CTEMP1-V2 (IPTR+J)  )  /CABS  (CTEflPl) 

IF  (CHANGE. GT.CHNMAX)  CHNM AX=CH ANGE 
CHNAVG=CHNAVG+CHANGE 
V2  (IPTR  +  J) =CT  EME 1 
V 1 (IPTR+J) =CZERC 
V3 (IPTR+J) =CZ  ERC 
250  CONTINUE 
260  CONTINUE 

ASK  WHETHER  OR  NOT  TO  STOP  AETER  REPORTING  *  FIELD  ERROR 

CHN A VG=  (CHNAVG+100. 0)  /  EE ACT 

CHNMAX=CHNMAX*  1C0.C 

WRITE  (3,1300) CHNAVG,CHNM AX 

R  FA D  (  1 ,  1100) FLAG  1 

CCONT=COUNT+1 

IF (FLAG1.EQ. 0)  GO  TC  105 

WRITE  (3,1400)  (Y  1  (I)  ,1=1,  NO) 

WRITE  (21,1700)  (Y1  (I)  ,1  =  1,  NO) 

WRITE  (3, 1400)  (Y2  (I)  ,1=1, NO) 

WRITE  (21,1  700)  (Y2  (I)  ,1  =  1,  NO) 

WRITE  (3,1400)  (Y3(I)  ,1=1, NO) 

WRITE  (21, 1700)  (Y  3  (I)  ,1=1,  NO) 

ASK  IF  FIELD  SHCULE  BE  PRINTED  OUT  ALSO 
WRITE  (3,1500) 

READ  (1,1100)  FLAG2 
IF  (FLAG2.NE. 0)  RETURN 
WRITE  (3,1600)  (72  (1)  ,1=1, M) 

RETURN 

1000  FORMAT (10E0.0) 

1100  FORMAT (410) 

1200  FORM  AT ( 1 H  ,'AVG  CURRENT  CHANGE=* , E 14.  7,  '  *•/ 

$  1 H  , ' MAX  CURRENT  CHANGF= » , E 1 4. 7, *  **// 

$  1 H  , 'AFTER*, 14, '  ITERATIONS*/) 
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1300  FORMAT  ( 1 H  ,'AVG  FIELD  ERROR  =  *,£15.7,*  %»/ 

$  1 H  ,  '  MAX  FIELD  ERROR  =  *,E15.7,*  *•/ 

$  1 H  , 'CONTINUE  ITERATIONS?  0  FOR  YES,  1  FOR  NO,  AND  RETURN*/) 

1  400  FORMAT  (1 H  , ' CDRRE NTS ' //  ( 1 H  ,10111-4)) 

1500  FORMAT  (1H  , 'PRINT  FIELDS?  0  FOR  YES,  1  FOR  NO,  THEN  RETURN'/) 
1600  FORMAT  (1H  , *  R  FSDLTA  NT  FIELDS '//( 1H  ,10E11.4)) 

1700  FORMAT  (8E  14.  6) 

END 

C 

C  COMPUTE  MUTUAL  IMPEDANCE  BETWEEN  TWO  PARALLEL  SEGMENTS 

C  OF  THIN  WIRE  DIPGIES  OF  SAME  LENGTH 

FUNCTION  ZMNG(Z1,Z2,LBY2,RSQ) 

REAL  LEY2 

COMPLEX  ZMNG,CMFLX 
DZ=AB  S  (Z1-Z2) 

CC=2. 0*CCS (LEY2) 

C5Q=CC*CC 
D  1=  DZ 

D2=DZ+LBY2 
D3=DZ-LBY2 
D  4=D2  *LB Y  2 
D5=D3-LBY2 

U 1  =  SQRT (RSQ  +  D1*D  1) +  D  1 
U2=SQRT  (RSQ+D2*D2)  +  D2 
U 3=SQRT  (RSQ+  D3*D3)  +D3 
U4=SQRT  {R SQ  +  D4*  E4)  +D4 
U5=SQHT  (RSQ+D5  +  D5) +D5 

V  1=RS0/U 1 

V  2=R  SQ/U  2 
V3=RSQ/U3 
V4=RSQ/U4 
V5=RSC/05 

CALL  SICI  (SU1,CU1,U1) 

CALL  SICI  (SU2,CU2,U2) 

CALL  SICI  (SU3,CU3,U3) 

CALL  SICI  (SU4,CU4,U4) 

CALL  SICI (SU5,C05,U5) 

CALL  SICI  (SV1,CV1,V1) 

CALL  SICI  (SV2,CV2,V2) 

CALL  SICI  <SV3,CV3,V3) 

CALL  SICI  (SV4,CV4,V4) 

CALL  SICI  {SV5, CVS, V5) 

S  1  =  SIN  (D  1 ) 

S  2=S I N  (D2) 

S3=SIN  (D3) 

S4=SIN  (D4 ) 

S5=SIN (D5) 

C1=COS  (D  1 ) 

C2=COS { D  2 ) 

C3=COS (D3) 

C4=COS (D4) 

C5=COS ( D  5 ) 

RL=  (2.0  +  CSQ)  *  (Cl*  (CU1+CV1)  +S 1  *  (SU  1-SV 1 )  ) 

$  -2. 0*CC*  (C2*  (CU2+CV2)  +  S2*  (SU2-SV2)  +C3*  (CU3+CV3) 

$  +  S3*  (SU3-3V3)  )  +C4*  (CU4+CV4)  +  S4*  (SU4-SV4) 

$  +  C5*  (C05+CV5) +S5*  (3U5-SV5) 

A  G=  (2.0+CSQ)*  (SI*  (CU1-CV1)-C1*(SU1*SV1)) 

$  -2.0*CC*(S2*(CU2-CV2)-C2*  (SU2+SV2) +S3* (C03-CV3) 

$  -C3*  (SD3  +  SV3) )  +S4*  (CU4-CV4) -C4*  (SU4  +  SV4) 

$  +S5*  (CU5-CV5)-C5*  (SU5+SV5) 
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ZMNG=15.0*CMPIX(Rl,AG)/  (SIN  (LBY2)  *SIN  (L9Y2)  ) 

RETURN 

END 

C 

C 

SUBROUTINE  SICI  ( SI ,CI , X) 

Z=ABS  (X) 

IF  (Z-4.  )  1,  1,4 
1  Y=  (4.-Z)*  (4.+Z) 

3  SI=X*(  <((  (  1. 753  141 E- 9*7*  1. 5  6  89  8  8  E— 7)  *Y  +  1.  3741  68E-5)  *T+6. 939889E-4) 
1*Y+1. 964882 E- 2) *  Y  +  4  . 395509E-1) 

CI={ (5. 772156E-1 +  AL0G (Z) ) /Z-Z* { ( ( (  ( 1 .386985E- 1 0*Y+ 1. 584996E-8) *Y 
1  +  1.  725752 E-6) *Y+1.  185999E-4) *Y+4. 99  0920E- 3) *Y+1. 315308E- 1)) *Z 
RETURN 

4  SI  =  SIN(Z) 

Y=COS  (Z) 

Z=4./Z 

U =<(({((( ( 4 . 048069E-3  *Z-2 . 279  143E-2) *Z+5. 515070E-2) *Z-7. 261642E-2) 
1*Z+4. 987716 E-2)*Z -3. 332519 E- 3)  *Z-2.  3  146 17E-2)  *2-1.  134958E-5)*Z 
2+ 6.250  01  IE- 2) *Z  +  2. 5 6 3989E-  1 0 
V=  <({{(((  (  (-5.  108699E-3  *Z+2 . 8  19  179  E-  2)  ♦Z-e.  537283E-2)  *Z 
1+7.902034E-2)  *Z-4 . 40  04 1  6E-2)  *Z-7. 945556E-3)  *Z  +  2.  60  1 293E- 2)  *Z 
2-3.  764  0  00E-4) *Z-3.  122  4 18 E- 2) *Z-6. 6 4644  IE-7) *Z  +  2. 500000E- 1 
CI  =  Z*  (SI*V-Y*U) 

SI=-Z*  (SI*U  +  Y*V)  +1.570796 

RETURN 
END 
C 

c 

c 

SUBROUTINE  VOLTU  (M2,  M3,  V) 

C  UNIFORM  EXCITATION  WITH  SPECIFIED  AMPLITUDE  AND  PHASE 

COMPLEX  V (1) ,CMPLX, VALUE 
M23=M2*M3 
READ  (5, 1)  AM,  PH 

1  FORMAT  (2F0.  0) 

R  A  C=PH*  3.  14159/180. 

VALUE=CMPLX  (AM*COS  (RAD)  ,  AM*SIN (F AD)  ) 

DO  2  1=1, M23 
V (I) =VALUE 

2  CONTINUE 
WRITE (5,3) AM, PH 

3  FORMAT  (///IX, 'UNIFORM  VOLTAGE  EXCITATION  -'//IX, 

S  *  MAGNITUDE  = • , 1 P 1 E20. 5, »  PHASE  = • , IP  1 E 20. 5/) 

RETURN 

END 

C 

C 

c 

SUBROUTINE  VOLTC  (M2 , M 3 , V) 

C  READ  IN  COMPLEX  NUMBERS  AS  VOLTAGE  FOR  EACH  DIPOLE 

COMPLEX  V(1) 

M23=M2*M3 

READ  (5,  1)  (V  (I)  ,1  =  1  ,  M23) 

1  FORMAT  (2F0.0) 

HR ITE {5, 2) 

2  FORMAT (///IX,' ARBITRARY  VOLTAGE  EXCITATION  -'//) 

RETURN 

END 

C 


41300 

41400 

41500 

41600 

41700 

41800 

41900 

42000 

42100 

42200 

42300 

42400 

42  500 

42600 

42  700 

42800 

42900 

43000 

43100 

43200 

433  00 

43400 

43500 

43600 

43700 

43800 

43  900 

44000 

44100 

44200 

44300 

44400 

44500 

44600 

44  700 

44800 

44900 

45000 

45100 

45200 

45300 

45400 

45500 

45600 

45700 

45800 

45900 

46000 

46100 

46200 

46300 

46400 

46500 

46600 

46700 

46800 

46900 

47000 

47100 

47200 


C 
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SUEROUTINE  STEER  (M2 , H 3 , DX , DZ r PP) 

C  PROGRESSIVE  PHASE  SHIFT  CN  EACH  DIPOLE  - 

C  STEERING  THE  MAIN  BEAM  IN  BCTH  DIFECTICN 

COMPLEX  PP  (1) ,CMPIX 
READ  (5,1)  HZ, EX 

1  FORMAT  (2F0.0) 

PI2=6. 2831853 
THX=  (RX)  +PI2/360. 

THZ=  (RZ) *P 1 2/3 60. 

L=0 

XK=-DX*COS (THX) 

ZK=-DZ*COS (THZ) 

DO  100  1=1, M3 
DO  100  J=1,M2 
L=L  + 1 

PH=FLO AT  (J-1)  *XK  +  FLOAT  (1-1)  *ZK 
PP  (L)  =  1.  0*CMP1X  (COS  (PH)  ,  SIN  (PH)  ) 

100  CONTINOE 

HRITE(5,20)  RX , RZ 

20  FORMAT  {///IX, ’BEAM  STEERING  =*,F10.5,'  DEGREES  IN  PHI  ANGLE', 
S/1X,14X,'=’,F10.5, *  DEGREES  IN  THE  ANGIE’//) 

RETORN 

ENE 

C 

C 

SUBROUTINE  T  AP  (M2  ,  M3  ,  EX ,  EZ ,  V V) 

C  MAGNITUDE  TAPER  OF  EXCITATION  IN  BOTH  DIRECTION 

COMPLEX  VV(1) 

PI2=6. 2831353 
HFX=DX*  (M2-  1)  *0.  5 
HFZ=DZ*  (M3-  1)  *0.5 
L=  0 

8R  7TE (5, 2) 

2  FORMAT  {///IX,  ’  EXPONENTIAL  TAPERED  IN  MAGNITUDE'///) 

DO  100  1=1, M3 

Z=  (  (1-1)  *DZ-HFZ)  / P 12 
FUNZ  =  EXP  (- ABS  (Z)  ) 

DO  100  J=  1  ,  M2 
L=1  +  1 

X=  (  (J—  1)*DX-HFX)/PI2 
VV  (L)  =EXP  (—ABS  (X)  )  *FUNZ 
100  CONTINUE 
RETURN 
ENE 
C 
C 
C 

SUEROUTINE  VOLTK  ( M 2 , M3 , V) 

COMPLEX  V  { 1 )  ,VK(20)  ,VJ{20) 

READ  (5,1)  (VK(I)  ,1=1, M3) 

READ  (5,1)  (VJ  (I)  ,1=1, M2) 

1  FORMAT (10F0.0) 

L-  0 

DO  2  1=1,  M3 
DO  2  J=  1 , M2 
L=I  + 1 

V  (L)  =VK  (I)  *VJ(J) 

2  CONTINUE 
HR ITE (5,4) 


47  300 

47400 

47500 

47600 

47700 

47800 

47900 

48000 

48  100 

48200 

48300 

48400 

48500 

48600 

48700 

48800 

48900 

49  COO 

49100 

49200 

49300 

49400 

49500 

49600 

49700 

49800 

49900 

50000 

50  100 

50200 

50300 

50400 

50500 

50600 

50700 

50800 

50900 

51  000 

51  100 

51200 

51300 

51400 

51500 

51600 

51700 

51800 

51900 

5200C 

52100 

52200 

52300 

52400 

52500 

52600 

52  700 

52800 

52900 

53000 

53100 

53200 
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4  FORMAT  {///IX,  •VOLTAGE  EXCITATION  -  SPECIFIED  BY  HOW, 
S'  AND  COLUMN' ,//) 

BETURN 

ENE 

C 

c 

SUFHOUTINE  PHASE (M2, M3, PP) 

C  PBOGBESSIV  E  PHASE  SHIFT  ON  EACH  EIPOLE  - 

COMPLEX  PP  (1)  ,CMPLX 
BEAD  (5, 1)  BZ,RX 
1  FORMAT  (2F0.  0) 

PI 2= 6. 2831853 
THX=  (BX)  +PI2/360. 

TH  Z—  (RZ) *P 12/360. 

L=  0 

DO  100  1=1, M3 
DO  100  J=  1  , M2 
L=I  + 1 

PH=FLO AT  (J-1) *THX  + FLOAT (1-1) *THZ 
PP  (L)  =  1.0+CMPLX  (COS  (PH)  ,STN  (PH)  ) 

100  CONTINUE 

WRITE  (5,20) RX,BZ 

20  FORM  AT  (///I  X, 'PROGRESSIVE  PHASE  SHIFT  =  »,F10.5, 
S'DEGREES  IN  ROW  DIE ECT ION '/25X , ' =  *,F10.5, 

S 'DEGREES  IN  COLUMN  DIRECTION  * ,///) 

RETURN 

ENE 

SUBROUTINE  SIZE  (LC ,  MC ,  LM,  MM  ,L ,  M  , N) 

L  =  LO*  3- 2 
M— MO*  3-2 
LTEMP-L 
M  T  E  M  P  =  M 
L  M=  0 
MM=0 

1  L=L/2 
Lfl=LM+  1 

IF  (L.GT. 1)  GO  TO  1 

2  M=M/2 
MM  =  MM  +  1 

IF  (M. GT. 1)  GO  TO  2 

L=2**1M 

M=2**MM 

IF (LTEMP.GT.L)  LN=IM+ 1 

IF  (LTEMP.GT.L)  L=L*2 

IF (MTFMP.GT. M)  MM  =  MM+  1 

IF  (MTEMP.GT.  M)  M=M*2 

N  =  M*L 

RETURN 

END 

SUBROUTINE  FFT  (X, M, START, ST EP) 

COMPLEX  X  ( 16384) ,U,K,T 
INTEGER  START, STEP, SEIFF 
N=2**N 

SEIFF=STEP-START 

N  V 2  =N/2*STEP 

N  M 1  =  (N-2) *STEP  +  STAFT 

N  =  (N-1)  +STEP  +  START 

J  =ST A  RT 

DO  8  I=START,NM1,STEP 

IF  (I-GE.  J)  GO  TO  5 
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53300 

T  =  X<J) 

53400 

X(J)  =  X(I) 

53500 

XII)  =T 

53600 

5 

K  =  NV  2 

53700 

6 

IF  (K-SDIFF.GE.  J)  GC  TC  7 

53800 

J  =  J-K 

53  900 

K  =K/2 

54000 

GC  TO  6 

54100 

7 

J  =  J  +  K 

54200 

8 

CONTINUE 

54300 

PI  =3.14159265358979 

54400 

DO  20  L=1,M 

54500 

LF  =  2**L 

54600 

LE 1=LE/2 

54700 

LSTEP=LE1*STEP 

54800 

U  ={1.0, 0.0) 

54  900 

A  NGLE=PI/FLC  AT  (LEI) 

55000 

W  =CMPLX  {CCS  (ANGLE)  ,-SIN  (ANGLE)  ) 

55100 

LEI  =LSTEP+START-STEP 

55200 

LE  =LE*STEP 

55300 

DO  20  J= START, LEI, STEP 

55400 

DC  10  1= J , N , I E 

55500 

IP  =I+LSTEP 

55600 

T  =X{IP)*U 

55700 

X  (IP)=X  (I)  -T 

55800 

X(I)  =X{I)+T 

55900 

10 

CONTINUE 

56000 

U=U*W 

56100 

20 

CONTINUE 

56200 

RETURN 

56300 

END 

56400 

C 

56500 

SUBROUTINE  TWODF {X, N, M, L , KM , LW) 

56600 

COMPLEX  X (  16384) 

56700 

INTEGER  START, STEP 

56800 

START=  1 

56900 

STEP  =1 

57000 

DC  10  1=1, M 

57100 

CALL  FFT (X,LM, START, STEF) 

57200 

START=START+L 

57300 

10 

CONTINUE 

57400 

STEP  =L 

57500 

DC  20  1=1, L 

57600 

ST  ART=I 

57700 

CALL  EFT (X, MM, START, STEP) 

57800 

20 

CONTINUE 

57  90C 

RETURN 

58000 

END 

58100 

SUBROUTINE  I  EFT  (X,  M  ,  ST  A  BT,  ST  EP) 

58200 

COMPLEX  X  (  16384) ,U,W,T 

58300 

INTEGER  START, STEP, SDIFF 

58400 

N  =  2**M 

58500 

S  CIFF=STEP-ST  ART 

58600 

N  V  2  =N/2*STEP 

58700 

N  M 1  =  (N-2) * ST EP  + START 

58800 

NEXP  =  (N-  1) +STEP  +  STABT 

58900 

J  =ST A  RT 

59000 

DO  8  I=START,NM1,STEP 

59100 

IE  (I.  GE.  J)  GC  TC  5 

59200 

T  =X  (J) 
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59300 

X(J)  =  X(I) 

59400 

X  fl)  =T 

59500 

5 

K  =  NV2 

59600 

6 

IF  {K-SDIFF.GE. J)  GC  TO  7 

59700 

J  =  J-K 

59800 

K  =  K/2 

59900 

GC  TO  6 

60000 

7 

J  =  J*K 

60100 

8 

CONTINUE 

60200 

PI  =3.14159265358979 

60300 

DO  20  L=1,M 

60400 

LI  =2**L 

60500 

LE  1  =  LE/2 

60600 

LSTEP=LE1*STEP 

60700 

U  =(1.0, 0.0) 

60800 

A  NGLE=PI/FLO AT  (LEI) 

60900 

H  =CMPLX  (COS  (ANGLE)  , SIN  (ANGIE)  ) 

61000 

LEI  =LSTEP+STABT-STEP 

61100 

LE  =LE* STEP 

61200 

DC  20  J=START, LEI, STEP 

61300 

DO  10  I=J, NEXP, LE 

61400 

IP  =1 +LSTEP 

61500 

T  =  X(IP)*U 

61600 

X  (IP)  =X  (I)  -T 

61700 

X  (I)  =  X  (I)  +T 

61800 

10 

CCNTINUE 

61900 

U  =  U*« 

62000 

20 

CONTINUE 

62100 

RFTURN 

62200 

END 

62300 

C 

62400 

SUBROUTINE  ITWODF (X, N, ,LH) 

62500 

COMPLEX  X  (  16364) 

62600 

INTEGER  START, STEP 

62700 

ST  ART=  1 

62800 

STEP  =1 

62900 

DO  10  1=1, M 

63000 

CALL  IFFT (X,LM, START, STEP) 

63100 

START=START+L 

63200 

10 

CONTINUE 

63300 

STEP  =L 

63400 

DO  20  1=1, L 

63500 

START  =1 

63600 

CALL  IFFT  (X, MM, START, STEP) 

63700 

20 

CONTINUE 

63800 

FN=FLGAT  (N) 

63900 

DO  30  1=1, N 

64  COO 

X  (I)  =  X  (I)/FN 

64  100 

30 

CONTINUE 

64200 

RETURN 

64300 

END 
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00100 

C 

00  200 

C 

ROUTINE  TO  COMPOTE  THE  FAR  FIELD  PATTERN 

00300 

SUBROUTINE  PATTRN  (LO,  MO,  NO ,  Y ,  DX  ,  DZ) 

00400 

INTEGER  LO, MO, NO, PTR 

00500 

REAL  DX,DZ, ALPHA, COSTHE, MAGE, DANG, T1,T2,T3,T4 

00600 

COMPLEX  Y  (1) ,AF,XPHASE,ZPHASE,DXPH,DZPH,CMPLX 

00700 

REAL  PAT (361) 

00800 

WRITE  (3,  1000) 

00900 

READ(1,  1100)  NPTS,  DANG,  IFLAG 

01000 

ALPHA=0. 0 

01  100 

DO  300  11=1, NPTS 

01200 

COSTHE=COS (ALPHA*.  1745329 3E-0 1 ) 

01300 

T3=DX*COST  HE 

01400 

T  1=COS  (T3) 

01500 

T2=SIN(T3) 

01550 

T4=-T2 

01600 

DXPH=CMPLX(T1,T2) 

01700 

D2PH=DXPH 

01  800 

AF=  (0.0,0.  0) 

01900 

ZPHASE={1. 0,0.0) 

01  950 

IF  (IFLAG.  NE.  0)  DZ  PH=CMPL X  (T 1  ,T4) 

02000 

DO  200  1=1, MO 

02100 

PTR= (1-1) *LO 

02200 

XPHASE=  (1. 0,0.  0) 

02300 

DO  100  J=  1 , LO 

02400 

A  F=  A  F  +  Y (PTR  +  J) *XPHASE*ZPHASE 

02500 

XPHASE=XPH ASE*DXPH 

02600 

100 

CONTINUE 

02700 

ZPRASE=ZPHASE*DZPH 

02800 

200 

CONTINUE 

02900 

MAGE=CAPS (AF) 

03000 

PAT  (II)  =M AGE 

03100 

WPITE  (3,  1200)  ALPHA, MAGE, AF 

03200 

ALPH A=ALPHA+DANG 

03300 

300 

CONTINUE 

03400 

WRITE  (2  1,  1  300)  (PAT  (II)  ,11=1, NPTS) 

03  500 

1000 

FORMAT  (1H  , * NPTS  AND  ANGLE  INCREMENT?*) 

03600 

1  100 

FORMAT (10, E0. 0,10) 

03700 

1200 

FORMAT  (1 H  ,F8.  1,3E12.4) 

03800 

1300 

FORMAT  (8E  1 1 . 4) 

03900 

RETURN 

04000 

END 

04100 

C 

04200 

c 

ROUTINE  TO  ZEROIZE  CORNERS 

04300 

SUBROUTINE  SHAPE (Y, NS, LO) 

04400 

COMPLEX  Y  ( 1)  ,CZERO 

04  5  00 

CZERO=(0. 0,0.0) 

04600 

DO  30  1=  1  ,  NS—  1 

04700 

IPTR=  (1-1)  *LO 

04  800 

DO  10  J= 1 , NS— I 

04900 

Y (IPTR+J) =CZERO 

05000 

10 

CONTINUE 

05100 

DO  20  J= 2* NS +1-1 , 3*NS-2 

05200 

Y (IPTR+J) =CZERO 

05300 

20 

CONTINUE 

05400 

30 

CONTINUE 

05  500 

DO  60  I=2*NS, 3*NS-2 

05600 

IPTB=  (1-1)  *LO 

05700 

DO  40  J= 1 , I-2*NS+ 1 

05800 

Y  (IPTR  +  J) =CZFRO 

05900 

06000 

06100 

06200 

06300 

06400 

06500 

06600 

06700 

06  800 

06  9  0  0 

07000 

07100 

07  200 

07300 

07400 

07500 

07600 

07700 

07800 

07900 

08  000 

08100 

08150 

08200 

08300 

08400 

08500 

08600 

08700 

08750 

08800 

08900 

09  000 

09100 

09200 

09300 

09400 

09500 

09600 

09700 

09  800 

09900 

10000 

10100 

10200 

10250 

10300 

10400 

10500 

10600 

10700 

10800 

10850 

10900 

10  925 

10950 

10975 

11000 

11100 


72 


40  CONTINOE 

DO  50  J=5*NS-T-2,3*NS-2 
Y  (IPTR  + J)  =CZERO 
50  CONTINOE 
60  CONTINOE 
RETORN 
END 
C 

C  ROOTINE  TO  COMPOTE  THE  FUR  FIELD  PATTERN (3  EXPANSIONS) 

SOBROOTINE  PATT3E (LO, MO , NO, Y, DX , DZ , WLBY2) 

INTEGER  LO,MO,NO,PTR,CPTR 

REAL  DX,DZ, ALPHA,C0STHE,MAGE,DANG,T1,T2,T3,T4 ,DEL 
COMPLEX  Y  (1) ,AF1,AF2, AF3, XPH ASE,ZP  H ASE,DXPH, DZPH ,CMPLI ,  AF 
COMPLEX  PHASE, FAC1,FAC3 
REAL  PAT  (361) 

WRITE (3, 1000) 

R  EAD (  1 , 1 100)  NPTS,DANG,IFLAG 
ALPHA=0. 0 
DO  300  11=1, NPTS 
COSTHE=COS (ALPHA*. 17453293E-01) 

T3=DX*COSTHE 
T1=C0S  (T3) 

T2=SIN  (T3) 

T4=-T2 

DXPH=CMPLX  (T 1  ,T2) 

DZPH=DXPH 
AF1=(0. 0,0.0) 

AF2=(0. 0,0.0) 

AF3= (0.0,0. 0) 

ZPHASE=(1. 0,0.0) 

IF  (IFLAG.  NE.  0)  DZPH=CNPLX  (T1  ,T4) 

DO  200  1=1, MO 
PTR=  (1-1)  *LO 
XPHASE=(1. 0,0.0) 

DO  100  J=1,LO 
P  HAS  E=XPHASE*ZPHASE 
CPTR=PTR+ J 

AF1=AF1+Y  (CPTR)  *PHASE 
CPTR=CPTR*NO 
AF2=AF2+Y  (CPTR)  *PHASE 
CPTR=CPTR+NO 
AF3=AF3+Y (CPTR)  *PHASE 
XPHAS  E=XPHASE*DXPH 
100  CONTINOE 

ZPHASE=ZPHASE*DZPH 
200  CONTINOE 

AF=AF1+AF2+AF3 
T3=  WLBY2/2.  0*COSTHE 
T1=C0S  (T3) 

T2=SIN  (T3) 

F AC1=CHPLX  (T 1  ,T2) 

T  2=-T2 

FAC3=CMPLX  (T 1  ,T2) 

IF (IFLAG. EQ. 0)  GO  TO  250 
AF=AF1*FAC1+AF2+AF3*FAC3 
IF(COSTHE.GE. 0.99999)  GO  TO  250 

AF=AF*  (T1-COS  (WLBY2/2. 0)  )  /  (COS  (WLBY2*COSTHE)  -COS  (WLBY2)  ) 
250  CONTINOE 

M AGE=C ABS (AF) 

PAT  (II)  =HAGE 


11200 

11300 

11400 

11500 

11600 

11700 

11800 

11900 

12000 

12100 

12200 

12300 

12400 

12450 

12500 

12600 

12605 

12610 

12620 

12630 

12640 

12650 

12660 

12670 

12720 

12760 

12800 

12820 

12900 

12930 

12960 

13000 

13100 

13130 

13160 

13200 

13300 

13400 

13500 

13600 

13700 

13800 

13900 

14000 

14100 

14200 

14300 

14400 
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WRITE  (3,1 200)  ALPHA, MAGE, AF 
ALPHA=ALPHA+DANG 
300  CONTINUE 

WRITE  (2 1,1  300)  (PAT  (II)  ,II=1,NPTS) 

1000  FORMAT ( 1 H  , *NPTS  AND  ANGLE  INCREMENT? ' ) 

1100  FORMAT (10, E0. 0,10) 

1200  FORMAT ( 1H  , F8. 1 , 3E 12. 4) 

1300  FORMAT  (8E 1 1. 4) 

RETURN 

END 

C 

C  MAIN  PROGRAM 

COMPLEX  V  (12300) 

INTEGER  CODE,FLAG,LO, MO, NO, SKIP 
OPEN (UNIT= 21, FILE=* PATTRN.DAT') 

OPEN (UNIT=22, FILE=’ MXCUR.DAT*) 

READ (1, 200)  LO, MO,NO,DX,DZ,WLBY2, SKIP 
10  IF  (SKIP.LE.  0)  GO  TO  20 

READ  (22,300)  (7  (I)  ,1=  1 ,  NO) 

READ  (22,  300)  (7(1)  ,1=1, NO) 

READ(22,300)  (7(1)  ,1=1, NO) 

SKIP=SKIP- 1 
GC  TO  10 
20  CONTINUE 

THOPI=6. 2831853 

DX=DX*TWOPI 

DZ=DZ*TWOPI 

WLBY2= WLB Y2+TWOPI 

READ  (22,  300)  (7  (I)  ,1=  1,  NO) 

READ  (22,300)  (7  (T)  ,I  =  NO*  1 ,  NO*  NO) 

R  EAD (22,300)  (7  (I) ,I=NO+NO+1,NO+NO+NO) 

CALL  PATT3E (LO, MO, NO, 7, DX, DZ , WLBY2) 

NSIDE=  (LO  +  2)  /3 

NCOLS=LO 

NROWS=MO 

R  FAD  (  1, 200)  CODE, FLAG 

IF  (CODE.  EQ.  1)  CALL  70LTU  (NCOLS,NROHS,  7) 

IF  (CODE. EQ. 2)  CALL  70LTC  (NCOLS , NROWS, 7) 

IF (CO  BE. EO. 3)  CALI  TAP (NCOLS , NROWS , DX, DZ , 7) 

IF  (CODE. EQ.  4)  CALL  70LTK ( NCOLS , NROWS, 7) 

IF  (CODE.  EQ.  5)  CALL  STEER  (NCOLS, NROWS, DX,DZ, 7) 
IF  (CODE.  EQ.  6)  CALL  PHASE  (NCOLS,  NROWS, 7) 

CALL  SHAPE  (7,NSIDE,L0) 

CALL  PATTRN (LO ,M0,N0,7,DX,DZ) 

STOP 

200  FORMAT (3I0,3F0.0,I0) 

300  FORMAT (8E14. 6) 

END 
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